18  Statistique univariée & Intervalles de confiance

On entend par statistique univariée l’étude d’une seule variable, que celle-ci soit continue (quantitative) ou catégorielle (qualitative). La statistique univariée fait partie de la statistique descriptive.

18.1 Exploration graphique

Une première approche consiste à explorer visuelle la variable d’intérêt, notamment à l’aide de l’interface proposée par esquisse (cf Section 17.4).

Nous indiquons ci-après le code correspond aux graphiques ggplot2 les plus courants.

library(ggplot2)

18.1.1 Variable continue

Un histogramme est la représentation graphique la plus commune pour représenter la distribution d’une variable, par exemple ici la longueur des pétales (variable Petal.Length) du fichier de données datasets::iris. Il s’obtient avec la géométrie ggplot2::geom_histogram().

ggplot(iris) +
  aes(x = Petal.Length) +
  geom_histogram()
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Figure 18.1: un histogramme simple
Astuce

Il faut noter qu’il nous a suffit d’associer simplement la variable Petal.Length à l’esthétique x, sans avoir eu besoin d’indiquer une variable pour l’esthétique y.

En fait, ggplot2 associe par défaut à toute géométrie une certaine statistique. Dans le cas de ggplot2::geom_histogram(), il s’agit de la statistique ggplot2::stat_bin() qui divise la variable continue en classes de même largeur et compte le nombre d’observation dans chacune. ggplot2::stat_bin() renvoie un certain nombre de variables calculées (la liste complète est indiquée dans la documentation dans la section Compute variables), dont la variable count qui correspond au nombre d’observations la classe. On peut associer cette variable calculée à une esthétique grâce à la fonction ggplot2::after_stat(), par exemple aes(y = after_stat(count)). Dans le cas présent, ce n’est pas nécessaire car ggplot2 fait cette association automatiquement si l’on n’a pas déjà attribué une variable à l’esthétique y.

On peut personnaliser la couleur de remplissage des rectangles en indiquant une valeur fixe pour l’esthétique fill dans l’appel de ggplot2::geom_histogram() (et non via la fonction ggplot2::aes() puisqu’il ne s’agit pas d’une variable du tableau de données). L’esthétique colour permet de spécifier la couleur du trait des rectangles. Enfin, le paramètre binwidth permet de spécifier la largeur des barres.

ggplot(iris) +
  aes(x = Petal.Length) +
  geom_histogram(
    fill ="lightblue", 
    colour = "black", 
    binwidth = 1
  ) +
  xlab("Longeur du pétale") +
  ylab("Effectifs")
Figure 18.2: un histogramme personnalisé

On peut alternativement indiquer un nombre de classes avec bins.

ggplot(iris) +
  aes(x = Petal.Length) +
  geom_histogram(bins = 10, colour = "black")
Figure 18.3: un histogramme en 10 classes

Une représentation alternative de la distribution d’une variable peut être obtenue avec une courbe de densité, dont la particularité est d’avoir une surface sous la courbe égale à 1. Une telle courbe s’obtient avec ggplot2::geom_density(). Le paramètre adjust permet d’ajuster le niveau de lissage de la courbe.

ggplot(iris) +
  aes(x = Petal.Length) +
  geom_density(adjust = .5)
Figure 18.4: une courbe de densité

18.1.2 Variable catégorielle

Pour représenter la répartition des effectifs parmi les modalités d’une variable catégorielle, on a souvent tendance à utiliser des diagrammes en secteurs (camemberts). Or, ce type de représentation graphique est très rarement appropriée : l’œil humain préfère comparer des longueurs plutôt que des surfaces1.

1 Voir en particulier https://www.data-to-viz.com/caveat/pie.html pour un exemple concret.

Dans certains contextes ou pour certaines présentations, on pourra éventuellement considérer un diagramme en donut, mais le plus souvent, rien ne vaut un bon vieux diagramme en barres avec ggplot2::geom_bar(). Prenons pour l’exemple la variable occup du jeu de données hdv2003 du package questionr.

data("hdv2003", package = "questionr")
ggplot(hdv2003) +
  aes(x = occup) +
  geom_bar()
Figure 18.5: un diagramme en barres simple
Astuce

Là encore, ggplot2 a calculé de lui-même le nombre d’observations de chaque modalité, en utilisant cette fois la statistique ggplot2::stat_count().

Si l’on souhaite représenter des pourcentages plutôt que des effectifs, le plus simple est d’avoir recours à la statistique ggstats::stat_prop() du package ggstats2. Pour appeler cette statistique, on utilisera simplement stat = "prop" dans les géométries concernées.

2 Cette statistique est également disponible via le package GGally.

Cette statistique, qui sera également bien utile pour des graphiques plus complexes, nécessite qu’on lui indique une esthétique by pour dans quels sous-groupes calculés des proportions. Ici, nous avons un seul groupe considéré et nous souhaitons des pourcentages du total. On indiquera simplement by = 1.

Pour formater l’axe vertical avec des pourcentages, on pourra avoir recours à la fonction scales::label_percent() que l’on appellera via ggplot2::scale_y_continuous().

library(ggstats)
ggplot(hdv2003) +
  aes(x = occup, y = after_stat(prop), by = 1) +
  geom_bar(stat = "prop") +
  scale_y_continuous(labels = scales::label_percent())
Figure 18.6: un diagramme en barres épuré

Pour une publication ou une communication, il ne faut surtout pas hésiter à épurer vos graphiques (less is better!), voire à trier les modalités en fonction de leur fréquence pour faciliter la lecture (ce qui se fait aisément avec forcats::fct_infreq()).

ggplot(hdv2003) +
  aes(x = forcats::fct_infreq(occup), 
      y = after_stat(prop), by = 1) +
  geom_bar(stat = "prop", 
           fill = "#4477AA", colour = "black") +
  geom_text(
    aes(label = after_stat(prop) |> 
          scales::percent(accuracy = .1)),
    stat = "prop",
    nudge_y = .02
  ) +
  theme_minimal() +
  theme(
    panel.grid = element_blank(),
    axis.text.y = element_blank()
  ) +
  xlab(NULL) + ylab(NULL) +
  ggtitle("Occupation des personnes enquêtées")
Figure 18.7: un diagramme en barres épuré

Le diaporama ci-dessous vous permet de visualiser chaque étape du code.

18.2 Tableaux et tris à plat

Le package gtsummary constitue l’une des boites à outils de l’analyste quantitativiste, car il permet de réaliser très facilement des tableaux quasiment publiables en l’état. En matière de statistique univariées, la fonction clé est gtsummary::tbl_summary().

Commençons avec un premier exemple rapide. On part d’un tableau de données et on indique, avec l’argument include, les variables à afficher dans le tableau statistique (si on n’indique rien, toutes les variables du tableau de données sont considérées). Il faut noter que l’argument include de gtsummary::tbl_summary() utilise la même syntaxe dite tidy select que dplyr::select() (cf. Section 8.2.1). On peut indiquer tout autant des variables catégorielles que des variables continues.

library(gtsummary)
hdv2003 |> 
  tbl_summary(include = c(age, occup))
Characteristic N = 2,0001
age 48 (35, 60)
occup
    Exerce une profession 1,049 (52%)
    Chomeur 134 (6.7%)
    Etudiant, eleve 94 (4.7%)
    Retraite 392 (20%)
    Retire des affaires 77 (3.9%)
    Au foyer 171 (8.6%)
    Autre inactif 83 (4.2%)
1 Median (IQR); n (%)
Table 18.1: un tableau simple
Remarque sur les types de variables et les sélecteurs associés

gtsummary permets de réaliser des tableaux statistiques combinant plusieurs variables, l’affichage des résultats pouvant dépendre du type de variables.

Par défaut, gtsummary considère qu’une variable est catégorielle s’il s’agit d’un facteur, d’une variable textuelle ou d’une variable numérique ayant moins de 10 valeurs différentes.

Une variable sera considérée comme dichotomique (variable catégorielle à seulement deux modalités) s’il s’agit d’un vecteur logique (TRUE/FALSE), d’une variable textuelle codée yes/no ou d’une variable numérique codée 0/1.

Dans les autres cas, une variable numérique sera considérée comme continue.

Si vous utilisez des vecteurs labellisés (cf. Chapitre 12), vous devez les convertir, en amont, en facteurs ou en variables numériques. Voir l’extension labelled et les fonctions labelled::to_factor(), labelled::unlabelled() et unclass().

Au besoin, il est possible de forcer le type d’une variable avec l’argument type de gtsummary::tbl_summary().

gtsummary fournit des sélecteurs qui peuvent être utilisés dans les options des différentes fonctions, en particulier gtsummary::all_continuous() pour les variables continues, gtsummary::all_dichotolous() pour les variables dichotomiques et gtsummary::all_categorical() pour les variables catégorielles. Cela inclue les variables dichotomiques. Il faut utiliser all_categorical(dichotomous = FALSE) pour sélectionner les variables catégorielles en excluant les variables dichotomiques.

18.2.1 Thème du tableau

gtsummary fournit plusieurs fonctions préfixées theme_gtsummary_*() permettant de modifier l’affichage par défaut des tableaux. Vous aurez notez que, par défaut, gtsummary est anglophone.

La fonction gtsummary::theme_gtsummary_journal() permets d’adopter les standards de certaines grandes revues scientifiques telles que JAMA (Journal of the American Medical Association), The Lancet ou encore le NEJM (New England Journal of Medicine).

La fonction gtsummary::theme_gtsummary_language() permet de modifier la langue utilisée par défaut dans les tableaux. Les options decimal.mark et big.mark permettent de définir respectivement le séparateur de décimales et le séparateur des milliers. Ainsi, pour présenter un tableau en français, on appliquera en début de script :

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

Ce thème sera appliqué à tous les tableaux ultérieurs.

hdv2003 |> 
  tbl_summary(include = c(age, occup))
Caractéristique N = 2 0001
age 48 (35 – 60)
occup
    Exerce une profession 1 049 (52%)
    Chomeur 134 (6,7%)
    Etudiant, eleve 94 (4,7%)
    Retraite 392 (20%)
    Retire des affaires 77 (3,9%)
    Au foyer 171 (8,6%)
    Autre inactif 83 (4,2%)
1 Médiane (EI); n (%)
Table 18.2: un tableau simple en français

18.2.2 Étiquettes des variables

gtsummary, par défaut, prends en compte les étiquettes de variables (cf. Chapitre 11), si elles existent, et sinon utilisera le nom de chaque variable dans le tableau. Pour rappel, les étiquettes de variables peuvent être manipulées avec l’extension labelled et les fonctions labelled::var_label() et labelled::set_variable_labels().

Il est aussi possible d’utiliser l’option label de gtsummary::tbl_summary() pour indiquer des étiquettes personnalisées.

hdv2003 |> 
  labelled::set_variable_labels(
    occup = "Occupation actuelle"
  ) |> 
  tbl_summary(
    include = c(age, occup, heures.tv),
    label = list(age ~ "Âge médian")
  )
Caractéristique N = 2 0001
Âge médian 48 (35 – 60)
Occupation actuelle
    Exerce une profession 1 049 (52%)
    Chomeur 134 (6,7%)
    Etudiant, eleve 94 (4,7%)
    Retraite 392 (20%)
    Retire des affaires 77 (3,9%)
    Au foyer 171 (8,6%)
    Autre inactif 83 (4,2%)
heures.tv 2,00 (1,00 – 3,00)
    Manquant 5
1 Médiane (EI); n (%)
Table 18.3: un tableau étiquetté

Pour modifier les modalités d’une variable catégorielle, il faut modifier en amont les niveaux du facteur correspondant.

Remarque sur la syntaxe des options

De nombreuses options des fonctions de gtsummary peuvent s’appliquer seulement à une ou certaines variables. Pour ces options-là, gtsummary attends une formule de la forme variables concernées ~ valeur de l'option ou bien une liste de formules ayant cette forme.

Par exemple, pour modifier l’étiquette associée à une certaine variable, on peut utiliser l’option label de gtsummary::tbl_summary().

trial |> 
  tbl_summary(label = age ~ "Âge")

Lorsque l’on souhaite passer plusieurs options pour plusieurs variables différentes, on utilisera une list().

trial |> 
  tbl_summary(label = list(age ~ "Âge", trt ~ "Traitement"))

gtsummary est très flexible sur la manière d’indiquer la ou les variables concernées. Il peut s’agir du nom de la variable, d’une chaîne de caractères contenant le nom de la variable, ou d’un vecteur contenant le nom de la variable. Les syntaxes ci-dessous sont ainsi équivalentes.

trial |> 
  tbl_summary(label = age ~ "Âge")
trial |> 
  tbl_summary(label = "age" ~ "Âge")
v <- "age"
trial |> 
  tbl_summary(label = v ~ "Âge")

Pour appliquer le même changement à plusieurs variables, plusieurs syntaxes sont acceptées pour lister plusieurs variables.

trial |> 
  tbl_summary(label = c("age", "trt") ~ "Une même étiquette")
trial |> 
  tbl_summary(label = c(age, trt) ~ "Une même étiquette")

Il est également possible d’utiliser la syntaxe tidyselect et les sélecteurs de tidyselect comme tidyselect::everything(), tidyselect::starts_with(), tidyselect::contains() ou tidyselect::all_of(). Ces différents sélecteurs peuvent être combinés au sein d’un c().

trial |> 
  tbl_summary(
    label = everything() ~ "Une même étiquette"
  )
trial |> 
  tbl_summary(
    label = starts_with("a") ~ "Une même étiquette"
  )
trial |> 
  tbl_summary(
    label = c(everything(), -age, -trt) ~ "Une même étiquette"
  )
trial |> 
  tbl_summary(
    label = age:trt ~ "Une même étiquette"
  )

Bien sûr, il est possible d’utiliser les sélecteurs propres à gtsummary.

trial |> 
  tbl_summary(
    label = all_continuous() ~ "Une même étiquette"
  )
trial |> 
  tbl_summary(
    label = list(
      all_continuous() ~ "Variable continue",
      all_dichotomous() ~ "Variable dichotomique",
      all_categorical(dichotomous = FALSE) ~ "Variable catégorielle"
    )
  )

Enfin, si l’on ne précise rien à gauche du ~, ce sera considéré comme équivalent à everything(). Les deux syntaxes ci-dessous sont donc équivalentes.

trial |> 
  tbl_summary(label = ~ "Une même étiquette")
trial |> 
  tbl_summary(
    label = everything() ~ "Une même étiquette"
  )

18.2.3 Statistiques affichées

Le paramètre statistic permets de sélectionner les statistiques à afficher pour chaque variable. On indiquera une chaîne de caractères dont les différentes statistiques seront indiquées entre accolades ({}).

Pour une variable continue, on pourra utiliser {median} pour la médiane, {mean} pour la moyenne, {sd} pour l’écart type, {var} pour la variance, {min} pour le minimum, {max} pour le maximum, ou encore {p##} (en remplaçant ## par un nombre entier entre 00 et 100) pour le percentile correspondant (par exemple p25 et p75 pour le premier et le troisième quartile). Utilisez gtsummary::all_continous() pour sélectionner toutes les variables continues.

hdv2003 |>
  tbl_summary(
    include = c(age, heures.tv),
    statistic = 
      all_continuous() ~ "Moy. : {mean} [min-max : {min} - {max}]"
  )
Caractéristique N = 2 0001
age Moy. : 48 [min-max : 18 - 97]
heures.tv Moy. : 2,25 [min-max : 0,00 - 12,00]
    Manquant 5
1 Moy. : Moyenne [min-max : Étendue]
Table 18.4: statistiques personnalisées pour une variable continue

Il est possible d’afficher des statistiques différentes pour chaque variable.

hdv2003 |>
  tbl_summary(
    include = c(age, heures.tv),
    statistic = list(
      age ~ "Méd. : {median} [{p25} - {p75}]",
      heures.tv ~ "Moy. : {mean} ({sd})"
    )
  )
Caractéristique N = 2 0001
age Méd. : 48 [35 - 60]
heures.tv Moy. : 2,25 (1,78)
    Manquant 5
1 Méd. : Médiane [EI]; Moy. : Moyenne (ET)
Table 18.5: statistiques personnalisées pour une variable continue (2)

Pour les variables continues, il est également possible d’indiquer le nom d’une fonction personnalisée qui prends un vecteur et renvoie une valeur résumée. Par exemple, pour afficher la moyenne des carrés :

moy_carres <- function(x) {
  mean(x^2, na.rm = TRUE)
}
hdv2003 |>
  tbl_summary(
    include = heures.tv,
    statistic = ~ "MC : {moy_carres}"
  )
Caractéristique N = 2 0001
heures.tv MC : 8,20
    Manquant 5
1 MC : moy_carres
Table 18.6: statiques personnalisées pour une variable continue (3)

Pour une variable catégorielle, les statistiques possibles sont {n} le nombre d’observations, {N} le nombre total d’observations, et {p} le pourcentage correspondant. Utilisez gtsummary::all_categorical() pour sélectionner toutes les variables catégorielles.

hdv2003 |>
  tbl_summary(
    include = occup,
    statistic = all_categorical() ~ "{p} % ({n}/{N})"
  )
Caractéristique N = 2 0001
occup
    Exerce une profession 52 % (1 049/2 000)
    Chomeur 6,7 % (134/2 000)
    Etudiant, eleve 4,7 % (94/2 000)
    Retraite 20 % (392/2 000)
    Retire des affaires 3,9 % (77/2 000)
    Au foyer 8,6 % (171/2 000)
    Autre inactif 4,2 % (83/2 000)
1 % % (n/N)
Table 18.7: statiques personnalisées pour une variable catégorielle

Il est possible, pour une variable catégorielle, de trier les modalités de la plus fréquente à la moins fréquente avec le paramètre sort.

hdv2003 |>
  tbl_summary(
    include = occup,
    sort = all_categorical() ~ "frequency"
  )
Caractéristique N = 2 0001
occup
    Exerce une profession 1 049 (52%)
    Retraite 392 (20%)
    Au foyer 171 (8,6%)
    Chomeur 134 (6,7%)
    Etudiant, eleve 94 (4,7%)
    Autre inactif 83 (4,2%)
    Retire des affaires 77 (3,9%)
1 n (%)
Table 18.8: variable catégorielle triée par fréquence

Pour toutes les variables (catégorielles et continues), les statistiques suivantes sont également disponibles :

  • {N_obs} le nombre total d’observations,
  • {N_miss} le nombre d’observations manquantes (NA),
  • {N_nonmiss} le nombre d’observations non manquantes,
  • {p_miss} le pourcentage d’observations manquantes (i.e. N_miss / N_obs) et
  • {p_nonmiss} le pourcentage d’observations non manquantes (i.e. N_nonmiss / N_obs).

18.2.4 Affichage du nom des statistiques

Lorsque l’on affiche de multiples statistiques, la liste des statistiques est regroupée dans une note de tableau qui peut vite devenir un peu confuse.

tbl <- hdv2003 |>
  tbl_summary(
    include = c(age, heures.tv, occup),
    statistic = list(
      age ~ "{mean} ({sd})",
      heures.tv ~ "{median} [{p25} - {p75}]"
    )
  )
tbl
Caractéristique N = 2 0001
age 48 (17)
heures.tv 2,00 [1,00 - 3,00]
    Manquant 5
occup
    Exerce une profession 1 049 (52%)
    Chomeur 134 (6,7%)
    Etudiant, eleve 94 (4,7%)
    Retraite 392 (20%)
    Retire des affaires 77 (3,9%)
    Au foyer 171 (8,6%)
    Autre inactif 83 (4,2%)
1 Moyenne (ET); Médiane [EI]; n (%)
Table 18.9: tableau par défaut

La fonction gtsummary::add_stat_label() permets d’indiquer le type de statistique à côté du nom des variables ou bien dans une colonne dédiée, plutôt qu’en note de tableau.

tbl |> 
  add_stat_label()
Caractéristique N = 2 000
age, Moyenne (ET) 48 (17)
heures.tv, Médiane [EI] 2,00 [1,00 - 3,00]
    Manquant 5
occup, n (%)
    Exerce une profession 1 049 (52%)
    Chomeur 134 (6,7%)
    Etudiant, eleve 94 (4,7%)
    Retraite 392 (20%)
    Retire des affaires 77 (3,9%)
    Au foyer 171 (8,6%)
    Autre inactif 83 (4,2%)
Table 18.10: ajout du nom des statistiques 
tbl |> 
  add_stat_label(location = "column")
Caractéristique Statistique N = 2 000
age Moyenne (ET) 48 (17)
heures.tv Médiane [EI] 2,00 [1,00 - 3,00]
    Manquant n 5
occup

    Exerce une profession n (%) 1 049 (52%)
    Chomeur n (%) 134 (6,7%)
    Etudiant, eleve n (%) 94 (4,7%)
    Retraite n (%) 392 (20%)
    Retire des affaires n (%) 77 (3,9%)
    Au foyer n (%) 171 (8,6%)
    Autre inactif n (%) 83 (4,2%)
Table 18.11: ajout du nom des statistiques dans une colonne séparée

18.2.5 Forcer le type de variable

Comme évoqué plus haut, gtsummary détermine automatiquement le type de chaque variable. Par défaut, la variable age du tableau de données trial est traitée comme variable continue, death comme dichotomique (seule la valeur 1 est affichée) et grade comme variable catégorielle.

trial |>
  tbl_summary(
    include = c(grade, age, death)
  )
Caractéristique N = 2001
Grade
    I 68 (34%)
    II 68 (34%)
    III 64 (32%)
Age 47 (38 – 57)
    Manquant 11
Patient Died 112 (56%)
1 n (%); Médiane (EI)
Table 18.12: types de variable par défaut

Il est cependant possible de forcer un certain type avec l’argument type. Précision : lorsque l’on force une variable en dichotomique, il faut indiquer avec value la valeur à afficher (les autres sont alors masquées).

trial |>
  tbl_summary(
    include = c(grade, death),
    type = list(
      grade ~ "dichotomous",
      death ~ "categorical"
    ),
    value = grade ~ "III",
    label = grade ~ "Grade III"
  )
Caractéristique N = 2001
Grade III 64 (32%)
Patient Died
    0 88 (44%)
    1 112 (56%)
1 n (%)
Table 18.13: types de variable personnalisés

Il ne faut pas oublier que, par défaut, gtsummary traite les variables quantitatives avec moins de 10 valeurs comme des variables catégorielles. Prenons un exemple :

trial$alea <- sample(1:4, size = nrow(trial), replace = TRUE)
#| label: tbl-types-defaut-alea
#| tbl-cap: traitement par défaut d'une variable numérique à 4 valeurs uniques
trial |>
  tbl_summary(
    include = alea
  )
Caractéristique N = 2001
alea
    1 44 (22%)
    2 51 (26%)
    3 52 (26%)
    4 53 (27%)
1 n (%)

On pourra forcer le traitement de cette variable comme continue.

trial |>
  tbl_summary(
    include = alea,
    type = alea ~ "continuous"
  )
Caractéristique N = 2001
alea 3 (2 – 4)
1 Médiane (EI)
Table 18.14: forcer le traitement continu d’une variable numérique à 4 valeurs uniques

18.2.6 Afficher des statistiques sur plusieurs lignes (variables continues)

Pour les variables continues, gtsummary a introduit un type de variable "continuous2", qui doit être attribué manuellement via type, et qui permets d’afficher plusieurs lignes de statistiques (en indiquant plusieurs chaînes de caractères dans statistic). À noter le sélecteur dédié gtsummary::all_continuous2().

hdv2003 |>
  tbl_summary(
    include = c(age, heures.tv),
    type = age ~ "continuous2",
    statistic = 
      all_continuous2() ~ c(
        "{median} ({p25} - {p75})", 
        "{mean} ({sd})",
        "{min} - {max}"
      )
  )
Caractéristique N = 2 0001
age
    Médiane (EI) 48 (35 - 60)
    Moyenne (ET) 48 (17)
    Étendue 18 - 97
heures.tv 2,00 (1,00 – 3,00)
    Manquant 5
1 Médiane (EI)
Table 18.15: des statistiques sur plusieurs lignes (variables continues)

18.2.7 Mise en forme des statistiques

L’argument digits permet de spécifier comment mettre en forme les différentes statistiques. Le plus simple est d’indiquer le nombre de décimales à afficher. Il est important de tenir compte que plusieurs statistiques peuvent être affichées pour une même variable. On peut alors indiquer une valeur différente pour chaque statistique.

hdv2003 |>
  tbl_summary(
    include = c(age, occup),
    digits = list(
      all_continuous() ~ 1,
      all_categorical() ~ c(0, 1)
    )
  )
Caractéristique N = 2 0001
age 48,0 (35,0 – 60,0)
occup
    Exerce une profession 1 049 (52,5%)
    Chomeur 134 (6,7%)
    Etudiant, eleve 94 (4,7%)
    Retraite 392 (19,6%)
    Retire des affaires 77 (3,9%)
    Au foyer 171 (8,6%)
    Autre inactif 83 (4,2%)
1 Médiane (EI); n (%)
Table 18.16: personnalisation du nombre de décimales

Au lieu d’un nombre de décimales, on peut indiquer plutôt une fonction à appliquer pour mettre en forme le résultat. Par exemple, gtsummary fournit les fonctions suivantes : gtsummary::style_number() pour les nombres de manière générale, gtsummary::style_percent() pour les pourcentages (les valeurs sont multipliées par 100, mais le symbole % n’est pas ajouté), gtsummary::style_pvalue() pour les p-valeurs, gtsummary::style_sigfig() qui n’affiche, par défaut, que deux chiffres significatifs, ou encore gtsummary::style_ratio() qui est une variante de gtsummary::``style_sigfig() pour les ratios (comme les odds ratios) que l’on compare à 1.

Il faut bien noter que ce qui est attendu par digits, c’est une fonction et non le résultat d’une fonction. On indiquera donc le nom de la fonction sans parenthèse, comme dans l’exemple ci-dessous (même si pas forcément pertinent ;-)).

hdv2003 |>
  tbl_summary(
    include = age,
    digits = 
      all_continuous() ~ c(style_percent, style_sigfig, style_ratio)
  )
Caractéristique N = 2 0001
age 4 800 (35 – 60,0)
1 Médiane (EI)
Table 18.17: personnalisation de la mise en forme des nombres

Comme digits s’attend à recevoir une fonction (et non le résultat) d’une fonction, on ne peut pas passer directement des arguments aux fonctions style_*() de gtsummary. Pour cela il faut créer une fonction à la levée :

trial |>
  tbl_summary(
    include = marker,
    statistic = ~ "{mean} pour 100",
    digits = ~ function(x){style_percent(x, digits = 1)}
  )
Caractéristique N = 2001
Marker Level (ng/mL) 91,6 pour 100
    Manquant 10
1 Moyenne pour 100
Table 18.18: passer une fonction personnalisée à digits (syntaxe 1)

Depuis R 4.1, il existe une syntaxe raccourcie équivalente, avec le symbole \ à la place de function.

trial |>
  tbl_summary(
    include = marker,
    statistic = ~ "{mean} pour 100",
    digits = ~ \(x){style_percent(x, digits = 1)}
  )
Caractéristique N = 2001
Marker Level (ng/mL) 91,6 pour 100
    Manquant 10
1 Moyenne pour 100
Table 18.19: passer une fonction personnalisée à digits (syntaxe 2)

Une syntaxe alternative consiste à avoir recours à la fonction purrr::partial() qui permet d’appeler partiellement une fonction et de renvoyer une nouvelle fonction.

trial |>
  tbl_summary(
    include = marker,
    statistic = ~ "{mean} pour 100",
    digits = ~ purrr::partial(style_percent, digits = 1)
  )
Caractéristique N = 2001
Marker Level (ng/mL) 91,6 pour 100
    Manquant 10
1 Moyenne pour 100
Table 18.20: passer une fonction personnalisée à digits (syntaxe 3)

À noter dans l’exemple précédent que les fonctions style_*() de gtsummary tiennent compte du thème défini (ici la virgule comme séparateur de décimale).

Pour une mise en forme plus avancée des nombres, il faut se tourner vers l’extension scales et ses diverses fonctions de mise en forme comme scales::label_number() ou scales::label_percent().

ATTENTION : les fonctions de scales n’héritent pas des paramètres du thème gtsummary actif. Il faut donc personnaliser le séparateur de décimal dans l’appel à la fonction.

trial |>
  tbl_summary(
    include = marker,
    statistic = ~ "{mean}",
    digits = ~ scales::label_number(
      accuracy = .01, 
      suffix = " ng/mL", 
      decimal.mark = ","
    )
  )
Caractéristique N = 2001
Marker Level (ng/mL) 0,92 ng/mL
    Manquant 10
1 Moyenne
Table 18.21: passer une fonction personnalisée à digits (syntaxe 4)

18.2.8 Données manquantes

Le paramètre missing permets d’indiquer s’il faut afficher le nombre d’observations manquantes (c’est-à-dire égales à NA) : "ifany" (valeur par défaut) affiche ce nombre seulement s’il y en a, "no" masque ce nombre et "always" force l’affichage de ce nombre même s’il n’y pas de valeur manquante. Le paramètre missing_text permets de personnaliser le texte affiché.

hdv2003 |>
  tbl_summary(
    include = c(age, heures.tv),
    missing = "always",
    missing_text = "Nbre observations manquantes"
  )
Caractéristique N = 2 0001
age 48 (35 – 60)
    Nbre observations manquantes 0
heures.tv 2,00 (1,00 – 3,00)
    Nbre observations manquantes 5
1 Médiane (EI)
Table 18.22: forcer l’affichage des valeurs manquantes

Il est à noter, pour les variables catégorielles, que les valeurs manquantes ne sont jamais pris en compte pour le calcul des pourcentages. Pour les inclure dans le calcul, il faut les transformer en valeurs explicites, par exemple avec forcats::fct_na_value_to_level() de forcats.

hdv2003 |>
  dplyr::mutate(
    trav.imp.explicit = trav.imp |> 
      forcats::fct_na_value_to_level("(non renseigné)")
  ) |> 
  tbl_summary(
    include = c(trav.imp, trav.imp.explicit)
  )
Caractéristique N = 2 0001
trav.imp
    Le plus important 29 (2,8%)
    Aussi important que le reste 259 (25%)
    Moins important que le reste 708 (68%)
    Peu important 52 (5,0%)
    Manquant 952
trav.imp.explicit
    Le plus important 29 (1,5%)
    Aussi important que le reste 259 (13%)
    Moins important que le reste 708 (35%)
    Peu important 52 (2,6%)
    (non renseigné) 952 (48%)
1 n (%)
Table 18.23: valeurs manquantes explicites (variable catégorielle)

18.2.9 Ajouter les effectifs observés

Lorsque l’on masque les manquants, il peut être pertinent d’ajouter une colonne avec les effectifs observés pour chaque variable à l’aide de la fonction gtsummary::add_n().

hdv2003 |>
  tbl_summary(
    include = c(heures.tv, trav.imp),
    missing = "no"
  ) |> 
  add_n()
Caractéristique N N = 2 0001
heures.tv 1 995 2,00 (1,00 – 3,00)
trav.imp 1 048
    Le plus important
29 (2,8%)
    Aussi important que le reste
259 (25%)
    Moins important que le reste
708 (68%)
    Peu important
52 (5,0%)
1 Médiane (EI); n (%)
Table 18.24: ajouter une colonne avec les effectifs observés

18.3 Calcul manuel

18.3.1 Variable continue

R fournit de base toutes les fonctions nécessaires pour le calcul des différentes statistiques descriptives :

Si la variable contient des valeurs manquantes (NA), ces fonctions renverront une valeur manquante, sauf si on leur précise na.rm = TRUE.

hdv2003$heures.tv |> mean()
[1] NA
hdv2003$heures.tv |> mean(na.rm = TRUE)
[1] 2.246566
hdv2003$heures.tv |> sd(na.rm = TRUE)
[1] 1.775853
hdv2003$heures.tv |> min(na.rm = TRUE)
[1] 0
hdv2003$heures.tv |> max(na.rm = TRUE)
[1] 12
hdv2003$heures.tv |> range(na.rm = TRUE)
[1]  0 12
hdv2003$heures.tv |> median(na.rm = TRUE)
[1] 2

La fonction quantile() permets de calculer tous types de quantiles.

hdv2003$heures.tv |> quantile(na.rm = TRUE)
  0%  25%  50%  75% 100% 
   0    1    2    3   12 
hdv2003$heures.tv |> 
  quantile(
    probs = c(.2, .4, .6, .8),
    na.rm = TRUE
  )
20% 40% 60% 80% 
  1   2   2   3

La fonction summary() renvoie la plupart de ces indicateurs en une seule fois, ainsi que le nombre de valeurs manquantes.

hdv2003$heures.tv |> summary()
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
  0.000   1.000   2.000   2.247   3.000  12.000       5

18.3.2 Variable catégorielle

Les fonctions de base pour le calcul d’un tri à plat sont les fonctions table() et xtabs(). Leur syntaxe est quelque peu différente. On passe un vecteur entier à table() alors que la syntaxe de xtabs() se rapproche de celle d’un modèle linéaire : on décrit le tableau attendu à l’aide d’une formule et on indique le tableau de données. Les deux fonctions renvoient le même résultat.

tbl <- hdv2003$trav.imp |> table()
tbl <- xtabs(~ trav.imp, data = hdv2003)
tbl <- hdv2003 |> xtabs(~ trav.imp, data = _)
tbl
trav.imp
           Le plus important Aussi important que le reste 
                          29                          259 
Moins important que le reste                Peu important 
                         708                           52

Comme on le voit, il s’agit du tableau brut des effectifs, sans les valeurs manquantes, et pas vraiment lisible dans la console de R.

Pour calculer les proportions, on appliquera prop.table() sur la table des effectifs bruts.

prop.table(tbl)
trav.imp
           Le plus important Aussi important que le reste 
                  0.02767176                   0.24713740 
Moins important que le reste                Peu important 
                  0.67557252                   0.04961832

Pour la réalisation rapide d’un tri à plat, on pourra donc préférer la fonction questionr::freq() qui affiche également le nombre de valeurs manquantes et les pourcentages, en un seul appel.

hdv2003$trav.imp |> 
  questionr::freq(total = TRUE)
                                n     %  val%
Le plus important              29   1.5   2.8
Aussi important que le reste  259  13.0  24.7
Moins important que le reste  708  35.4  67.6
Peu important                  52   2.6   5.0
NA                            952  47.6    NA
Total                        2000 100.0 100.0

18.4 Intervalles de confiance

18.4.1 Avec gtsummary

La fonction gtsummary::add_ci() permet d’ajouter des intervalles de confiance à un tableau créé avec gtsummary::tbl_summary().

Avertissement

Par défaut, pour les variables continues, gtsummary::tbl_summary() affiche la médiane tandis que gtsummary::add_ci() calcule l’intervalle de confiance d’une moyenne !

Il faut donc :

  • soit afficher la moyenne dans gtsummary::tbl_summary() à l’aide du paramètre statistic ;
  • soit calculer les intervalles de confiance d’une médiane (méthode "wilcox.text") via le paramètre method de gtsummary::add_ci().
hdv2003 |>
  tbl_summary(
    include = c(age, heures.tv, trav.imp),
    statistic = age ~ "{mean} ({sd})"
  ) |> 
  add_ci(
    method = heures.tv ~ "wilcox.test"
  )
Caractéristique N = 2 0001 95% CI2
age 48 (17) 47, 49
heures.tv 2,00 (1,00 – 3,00) 2,5, 2,5
    Manquant 5
trav.imp

    Le plus important 29 (2,8%) 1,9%, 4,0%
    Aussi important que le reste 259 (25%) 22%, 27%
    Moins important que le reste 708 (68%) 65%, 70%
    Peu important 52 (5,0%) 3,8%, 6,5%
    Manquant 952
1 Moyenne (ET); Médiane (EI); n (%)
2 IC = intervalle de confiance
Table 18.25: ajouter les intervalles de confiance

L’argument statistic permet de personnaliser la présentation de l’intervalle ; conf.level de changer le niveau de confiance et style_fun de modifier la mise en forme des nombres de l’intervalle.

hdv2003 |>
  tbl_summary(
    include = c(age, heures.tv),
    statistic = ~ "{mean}"
  ) |> 
  add_ci(
    statistic = ~ "entre {conf.low} et {conf.high}",
    conf.level = .9,
    style_fun = ~ purrr::partial(style_number, digits = 1)
  )
Caractéristique N = 2 0001 90% CI2
age 48 entre 47,5 et 48,8
heures.tv 2,25 entre 2,2 et 2,3
    Manquant 5
1 Moyenne
2 IC = intervalle de confiance
Table 18.26: des intervalles de confiance personnalisés

18.4.2 Calcul manuel

Le calcul de l’intervalle de confiance d’une moyenne s’effectue avec la fonction t.test().

hdv2003$age |> t.test()

    One Sample t-test

data:  hdv2003$age
t = 127.12, df = 1999, p-value < 2.2e-16
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
 47.41406 48.89994
sample estimates:
mean of x 
   48.157

Le résultat renvoyé est une liste contenant de multiples informations.

hdv2003$age |> t.test() |> str()
List of 10
 $ statistic  : Named num 127
  ..- attr(*, "names")= chr "t"
 $ parameter  : Named num 1999
  ..- attr(*, "names")= chr "df"
 $ p.value    : num 0
 $ conf.int   : num [1:2] 47.4 48.9
  ..- attr(*, "conf.level")= num 0.95
 $ estimate   : Named num 48.2
  ..- attr(*, "names")= chr "mean of x"
 $ null.value : Named num 0
  ..- attr(*, "names")= chr "mean"
 $ stderr     : num 0.379
 $ alternative: chr "two.sided"
 $ method     : chr "One Sample t-test"
 $ data.name  : chr "hdv2003$age"
 - attr(*, "class")= chr "htest"

Si l’on a besoin d’accéder spécifiquement à l’intervalle de confiance calculé :

hdv2003$age |> t.test() |> purrr::pluck("conf.int")
[1] 47.41406 48.89994
attr(,"conf.level")
[1] 0.95

Pour celui d’une médiane, on utilisera wilcox.test() en précisant conf.int = TRUE.

hdv2003$age |> wilcox.test(conf.int = TRUE)

    Wilcoxon signed rank test with continuity correction

data:  hdv2003$age
V = 2001000, p-value < 2.2e-16
alternative hypothesis: true location is not equal to 0
95 percent confidence interval:
 47.00001 48.50007
sample estimates:
(pseudo)median 
      47.99996 
hdv2003$age |> 
  wilcox.test(conf.int = TRUE) |>
  purrr::pluck("conf.int")
[1] 47.00001 48.50007
attr(,"conf.level")
[1] 0.95

Pour une proportion, on utilisera prop.test() en lui transmettant le nombre de succès et le nombre d’observations, qu’il faudra donc avoir calculé au préalable. On peut également passer une table à deux entrées avec le nombre de succès puis le nombre d’échecs.

Ainsi, pour obtenir l’intervalle de confiance de la proportion des enquêtés qui considèrent leur travail comme peu important, en tenant compte des valeurs manquantes, le plus simple est d’effectuer le code suivant3 :

3 Notez l’utilisation de rev() pour inverser le tableau créé avec xtabs() afin que le nombre de succès (TRUE) soit indiqués avant le nombre d’échecs (FALSE).

xtabs(~ I(hdv2003$trav.imp == "Peu important"), data = hdv2003) |> 
  rev() |> 
  prop.test()

    1-sample proportions test with continuity correction

data:  rev(xtabs(~I(hdv2003$trav.imp == "Peu important"), data = hdv2003)), null probability 0.5
X-squared = 848.52, df = 1, p-value < 2.2e-16
alternative hypothesis: true p is not equal to 0.5
95 percent confidence interval:
 0.03762112 0.06502346
sample estimates:
         p 
0.04961832

Par défaut, prop.test() produit un intervalle de confiance bilatéral en utilisant la méthode de Wilson avec correction de continuité. Pour plus d’information sur les différentes manières de calculer l’intervalle de confiance d’une proportion, on pourra se référer à ce billet de blog.

Astuce

Comme on le voit, il n’est pas aisé, avec les fonctions de R base de calculer les intervalles de confiance pour toutes les modalités d’une variable catégorielle.

On pourra éventuellement avoir recours à la petite fonction suivante qui réalise le tri à plat d’une variable catégorielle, calcule les proportions et leurs intervalles de confiance.

prop_ci <- function(x, conf.level = .95, correct = TRUE) {
  tbl <- as.data.frame(table(x), responseName = "n")
  tbl$N <- sum(tbl$n)
  tbl$prop <- tbl$n / tbl$N
  tbl$conf.low <- NA_real_
  tbl$conf.high <- NA_real_
  for (i in 1:nrow(tbl)) {
    test <- prop.test(
      x = tbl$n[i],
      n = tbl$N[i],
      conf.level = conf.level,
      correct = correct
    )
    tbl$conf.low[i] <- test$conf.int[1]
    tbl$conf.high[i] <- test$conf.int[2]
  }
  tbl
}
prop_ci(hdv2003$trav.imp)
                             x   n    N       prop   conf.low  conf.high
1            Le plus important  29 1048 0.02767176 0.01894147 0.04001505
2 Aussi important que le reste 259 1048 0.24713740 0.22151849 0.27463695
3 Moins important que le reste 708 1048 0.67557252 0.64614566 0.70369541
4                Peu important  52 1048 0.04961832 0.03762112 0.06502346

18.5 webin-R

La statistique univariée est présentée dans le webin-R #03 (statistiques descriptives avec gtsummary et esquisse) sur YouTube.