19  Statistique bivariée & Tests de comparaison

19.1 Deux variables catégorielles

19.1.1 Tableau croisé avec gtsummary

Pour regarder le lien entre deux variables catégorielles, l’approche la plus fréquente consiste à réaliser un tableau croisé, ce qui s’obtient très facilement avec l’argument by de la fonction gtsummary::tbl_summary() que nous avons déjà abordée dans le chapitre sur la statistique univariée (cf. Section 18.2).

Prenons pour exemple le jeu de données gtsummary::trial et croisons les variables stage et grade. On indique à by la variable à représenter en colonnes et à include celle à représenter en lignes.

library(gtsummary)
theme_gtsummary_language("fr", decimal.mark = ',')

Setting theme "language: fr"

trial |> 
  tbl_summary(
    include = stage,
    by = grade
  )

Caractéristique	I N = 68 1	II N = 68 1	III N = 64 1
T Stage
T1	17 (25%)	23 (34%)	13 (20%)
T2	18 (26%)	17 (25%)	19 (30%)
T3	18 (26%)	11 (16%)	14 (22%)
T4	15 (22%)	17 (25%)	18 (28%)
1 n (%)

Table 19.1: un tableau croisé avec des pourcentages en colonne

Par défaut, les pourcentages affichés correspondent à des pourcentages en colonne. On peut demander des pourcentages en ligne avec percent = "row" ou des pourcentages du total avec percent = "cell".

Il est possible de passer plusieurs variables à include mais une seule variable peut être transmise à by. La fonction gtsummary::add_overall() permet d’ajouter une colonne totale. Comme pour un tri à plat, on peut personnaliser les statistiques affichées avec statistic.

library(gtsummary)
trial |> 
  tbl_summary(
    include = c(stage, trt),
    by = grade,
    statistic = ~ "{p}% ({n}/{N})",
    percent = "row"
  ) |> 
  add_overall(last = TRUE)

Caractéristique	I N = 68 1	II N = 68 1	III N = 64 1	Overall N = 200 1
T Stage
T1	32% (17/53)	43% (23/53)	25% (13/53)	100% (53/53)
T2	33% (18/54)	31% (17/54)	35% (19/54)	100% (54/54)
T3	42% (18/43)	26% (11/43)	33% (14/43)	100% (43/43)
T4	30% (15/50)	34% (17/50)	36% (18/50)	100% (50/50)
Chemotherapy Treatment
Drug A	36% (35/98)	33% (32/98)	32% (31/98)	100% (98/98)
Drug B	32% (33/102)	35% (36/102)	32% (33/102)	100% (102/102)
1

Table 19.2: un tableau croisé avec des pourcentages en ligne

Important

Choisissez bien votre type de pourcentages (en lignes ou en colonnes). Si d’un point de vue purement statistique, ils permettent tous deux de décrire la relation entre les deux variables, ils ne correspondent au même story telling. Tout dépend donc du message que vous souhaitez faire passer, de l’histoire que vous souhaitez raconter.

gtsummary::tbl_summary() est bien adaptée dans le cadre d’une analyse de facteurs afin de représenter un outcome donné avec by et une liste de facteurs avec include.

Lorsque l’on ne croise que deux variables et que l’on souhaite un affichage un peu plus traditionnel d’un tableau croisé, on peut utiliser gtsummary::tbl_cross() à laquelle on transmettra une et une seule variable à row et une et une seule variable à col. Pour afficher des pourcentages, il faudra indiquer le type de pourcentages voulus avec percent.

trial |> 
  tbl_cross(
    row = stage,
    col = grade,
    percent = "row"
  )

	Grade			Total
	I	II	III
T Stage
T1	17 (32%)	23 (43%)	13 (25%)	53 (100%)
T2	18 (33%)	17 (31%)	19 (35%)	54 (100%)
T3	18 (42%)	11 (26%)	14 (33%)	43 (100%)
T4	15 (30%)	17 (34%)	18 (36%)	50 (100%)
Total	68 (34%)	68 (34%)	64 (32%)	200 (100%)

Table 19.3: un tableau croisé avec tbl_cross()

19.1.2 Représentations graphiques

La représentation graphique la plus commune pour le croisement de deux variables catégorielles est le diagramme en barres, que l’on réalise avec la géométrie ggplot2::geom_bar() et en utilisant les esthétiques x et fill pour représenter les deux variables.

library(ggplot2)
ggplot(trial) +
  aes(x = stage, fill = grade) +
  geom_bar() +
  labs(x = "T Stage", fill = "Grade", y = "Effectifs")

Figure 19.1: un graphique en barres croisant deux variables

On peut modifier la position des barres avec le paramètre position.

library(ggplot2)
ggplot(trial) +
  aes(x = stage, fill = grade) +
  geom_bar(position = "dodge") +
  labs(x = "T Stage", fill = "Grade", y = "Effectifs")

Figure 19.2: un graphique avec des barres côte à côte

Pour des barres cumulées, on aura recours à position = "fill". Pour que les étiquettes de l’axe des y soient représentées sous forme de pourcentages (i.e. 25% au lieu de 0.25), on aura recours à la fonction scales::percent() qui sera transmise à ggplot2::scale_y_continuous().

library(ggplot2)
ggplot(trial) +
  aes(x = stage, fill = grade) +
  geom_bar(position = "fill") +
  labs(x = "T Stage", fill = "Grade", y = "Proportion") +
  scale_y_continuous(labels = scales::percent)

Figure 19.3: un graphique en barres cumulées

Ajouter des étiquettes sur un diagramme en barres

Il est facile d’ajouter des étiquettes en ayant recours à ggplot2::geom_text(), à condition de lui passer les bons paramètres.

Tout d’abord, il faudra préciser stat = "count" pour indiquer que l’on souhaite utiliser la statistique ggplot2::stat_count() qui est celle utilisé par défaut par ggplot2::geom_bar(). C’est elle qui permets de compter le nombre d’observations.

Il faut ensuite utiliser l’esthétique label pour indiquer ce que l’on souhaite afficher comme étiquettes. La fonction after_stat(count) permet d’accéder à la variable count calculée par ggplot2::stat_count().

Enfin, il faut indiquer la position verticale avec ggplot2::position_stack(). En précisant un ajustement de vertical de 0.5, on indique que l’on souhaite positionner l’étiquette au milieu.

ggplot(trial) +
  aes(
    x = stage, fill = grade, 
    label = after_stat(count)
  ) +
  geom_bar() +
  geom_text(
    stat = "count", 
    position = position_stack(.5)
  )

Pour un graphique en barres cumulées, on peut utiliser de manière similaire ggplot2::position_fill(). On ne peut afficher directement les proportions avec ggplot2::stat_count(). Cependant, nous pouvons avoir recours à ggstats::stat_prop(), déjà évoquée dans le chapitre sur la statistique univariée (cf. Section 18.1.2) et dont le dénominateur doit être précisé via l’esthétique by.

library(ggstats)
ggplot(trial) +
  aes(
    x = stage, 
    fill = grade, 
    by = stage,
    label = scales::percent(after_stat(prop), accuracy = .1)
  ) +
  geom_bar(position = "fill") +
  geom_text(
    stat = "prop", 
    position = position_fill(.5)
  ) +
  scale_y_continuous(labels = scales::percent)

On peut aussi comparer facilement deux distributions, ici la proportion de chaque niveau de qualification au sein chaque sexe.

p <- ggplot(trial) +
  aes(
    x = stage,
    y = after_stat(prop),
    fill = grade, 
    by = grade,
    label = scales::percent(after_stat(prop), accuracy = 1)
  ) +
  geom_bar(
    stat = "prop", 
    position = position_dodge(.9)
  ) +
  geom_text(
    aes(y = after_stat(prop) - 0.01),
    stat = "prop", 
    position = position_dodge(.9),
    vjust = "top"
  ) +
  scale_y_continuous(labels = scales::percent)
p

Il est possible d’alléger le graphique en retirant des éléments superflus.

p + 
  theme_light() +
  xlab("") +
  ylab("") +
  labs(fill = "") +
  ggtitle("Distribution selon le niveau, par grade") +
  theme(
    panel.grid = element_blank(),
    panel.border = element_blank(),
    axis.text.y = element_blank(),
    axis.ticks = element_blank(),
    legend.position = "top"
  ) +
  scale_fill_brewer()

Le diaporama ci-dessous vous permet de visualiser chaque étape du code correspondant au graphique précédent.

19.1.3 Calcul manuel

Les deux fonctions de base permettant le calcul d’un tri à plat sont table() et xtabs() (cf. Section 18.3.2). Ces mêmes fonctions permettent le calcul du tri croisé de deux variables (ou plus). Pour table(), on passera les deux vecteurs à croisés, tandis que pour xtabs() on décrira le tableau attendu à l’aide d’une formule.

table(trial$stage, trial$grade)

    
      I II III
  T1 17 23  13
  T2 18 17  19
  T3 18 11  14
  T4 15 17  18

tab <- xtabs(~ stage + grade, data = trial)
tab

     grade
stage  I II III
   T1 17 23  13
   T2 18 17  19
   T3 18 11  14
   T4 15 17  18

Le tableau obtenu est basique et ne contient que les effectifs. La fonction addmargins() permet d’ajouter les totaux par ligne et par colonne.

tab |> addmargins()

     grade
stage   I  II III Sum
  T1   17  23  13  53
  T2   18  17  19  54
  T3   18  11  14  43
  T4   15  17  18  50
  Sum  68  68  64 200

Pour le calcul des pourcentages, le plus simple est d’avoir recours au package questionr qui fournit les fonctions questionr::cprop(), questionr::rprop() et questionr::prop() qui permettent de calculer, respectivement, les pourcentages en colonne, en ligne et totaux.

questionr::cprop(tab)

       grade
stage   I     II    III   Ensemble
  T1     25.0  33.8  20.3  26.5   
  T2     26.5  25.0  29.7  27.0   
  T3     26.5  16.2  21.9  21.5   
  T4     22.1  25.0  28.1  25.0   
  Total 100.0 100.0 100.0 100.0   

questionr::rprop(tab)

          grade
stage      I     II    III   Total
  T1        32.1  43.4  24.5 100.0
  T2        33.3  31.5  35.2 100.0
  T3        41.9  25.6  32.6 100.0
  T4        30.0  34.0  36.0 100.0
  Ensemble  34.0  34.0  32.0 100.0

questionr::prop(tab)

       grade
stage   I     II    III   Total
  T1      8.5  11.5   6.5  26.5
  T2      9.0   8.5   9.5  27.0
  T3      9.0   5.5   7.0  21.5
  T4      7.5   8.5   9.0  25.0
  Total  34.0  34.0  32.0 100.0

19.1.4 Test du Chi² et dérivés

Dans le cadre d’un tableau croisé, on peut tester l’existence d’un lien entre les modalités de deux variables, avec le très classique test du Chi² (parfois écrit χ² ou Chi²). Pour une présentation plus détaillée du test, on pourra se référer à ce cours de Julien Barnier.

Le test du Chi² peut se calculer très facilement avec la fonction chisq.test() appliquée au tableau obtenu avec table() ou xtabs().

tab <- xtabs(~ stage + grade, data = trial)
tab

     grade
stage  I II III
   T1 17 23  13
   T2 18 17  19
   T3 18 11  14
   T4 15 17  18

chisq.test(tab)

    Pearson's Chi-squared test

data:  tab
X-squared = 4.8049, df = 6, p-value = 0.5691

Si l’on est adepte de gtsummary, il suffit d’appliquer gtsummary::add_p() au tableau produit avec gtsummary::tbl_summary().

trial |> 
  tbl_summary(
    include = stage,
    by = grade
  ) |> 
  add_p()

Caractéristique	I N = 68 1	II N = 68 1	III N = 64 1	p-valeur 2
T Stage				0,6
T1	17 (25%)	23 (34%)	13 (20%)
T2	18 (26%)	17 (25%)	19 (30%)
T3	18 (26%)	11 (16%)	14 (22%)
T4	15 (22%)	17 (25%)	18 (28%)
1 n (%)
2 test du khi-deux d’indépendance

Table 19.4: un tableau croisé avec test du khi²

Dans notre exemple, les deux variables stage et grade ne sont clairement pas corrélées.

Un test alternatif est le test exact de Fisher. Il s’obtient aisément avec fisher.test() ou bien en le spécifiant via l’argument test de gtsummary::add_p().

tab <- xtabs(~ stage + grade, data = trial)
fisher.test(tab)

    Fisher's Exact Test for Count Data

data:  tab
p-value = 0.5801
alternative hypothesis: two.sided

trial |> 
  tbl_summary(
    include = stage,
    by = grade
  ) |> 
  add_p(test = all_categorical() ~ "fisher.test")

Caractéristique	I N = 68 1	II N = 68 1	III N = 64 1	p-valeur 2
T Stage				0,6
T1	17 (25%)	23 (34%)	13 (20%)
T2	18 (26%)	17 (25%)	19 (30%)
T3	18 (26%)	11 (16%)	14 (22%)
T4	15 (22%)	17 (25%)	18 (28%)
1 n (%)
2 test exact de Fisher

Table 19.5: un tableau croisé avec test exact de Fisher

Note

Formellement, le test de Fisher suppose que les marges du tableau (totaux lignes et colonnes) sont fixées, puisqu’il repose sur une loi hypergéométrique, et donc celui-ci se prête plus au cas des situations expérimentales (plans d’expérience, essais cliniques) qu’au cas des données tirées d’études observationnelles.

En pratique, le test du Chi² étant assez robuste quant aux déviations par rapport aux hypothèses d’applications du test (effectifs théoriques supérieurs ou égaux à 5), le test de Fisher présente en général peu d’intérêt dans le cas de l’analyse des tableaux de contingence.

19.1.5 Comparaison de deux proportions

Pour comparer deux proportions, la fonction de base est prop.test() à laquelle on passera un tableau à 2×2 dimensions.

tab <- xtabs(~ I(stage == "T1") + trt, data = trial)
tab |> questionr::cprop()

                trt
I(stage == "T1") Drug A Drug B Ensemble
           FALSE  71.4   75.5   73.5   
           TRUE   28.6   24.5   26.5   
           Total 100.0  100.0  100.0   

tab |> prop.test()

    2-sample test for equality of proportions with continuity correction

data:  tab
X-squared = 0.24047, df = 1, p-value = 0.6239
alternative hypothesis: two.sided
95 percent confidence interval:
 -0.2217278  0.1175050
sample estimates:
   prop 1    prop 2 
0.4761905 0.5283019 

Il est également envisageable d’avoir recours à un test exact de Fisher. Dans le cas d’un tableau à 2×2 dimensions, le test exact de Fisher ne teste pas si les deux proportions sont différents, mais plutôt si leur odds ratio (qui est d’ailleurs renvoyé par la fonction) est différent de 1.

fisher.test(tab)

    Fisher's Exact Test for Count Data

data:  tab
p-value = 0.5263
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
 0.4115109 1.5973635
sample estimates:
odds ratio 
 0.8125409 

Mais le plus simple reste encore d’avoir recours à gtsummary et à sa fonction gtsummary::add_difference() que l’on peut appliquer à un tableau où le paramètre by n’a que deux modalités. Pour la différence de proportions, il faut que les variables transmises à include soit dichotomiques.

trial |> 
  tbl_summary(
    by = trt,
    include = response
  ) |> 
  add_difference()

Caractéristique	Drug A N = 98 1	Drug B N = 102 1	Difference 2	95% IC 2,3	p-valeur 2
Tumor Response	28 (29%)	33 (34%)	-4,2%	-18% – 9,9%	0,6
Manquant	3	4
1 n (%)
2 2-sample test for equality of proportions with continuity correction
3 IC = intervalle de confiance

Table 19.6: différence entre deux proportions

Attention : si l’on passe une variable catégorielle à trois modalités ou plus, c’est la différence des moyennes standardisées (globale pour la variable) qui sera calculée et non la différence des proportions dans chaque groupe.

trial |> 
  tbl_summary(
    by = trt,
    include = grade
  ) |> 
  add_difference()

Caractéristique	Drug A N = 98 1	Drug B N = 102 1	Difference 2	95% IC 2,3
Grade			0,07	-0,20 – 0,35
I	35 (36%)	33 (32%)
II	32 (33%)	36 (35%)
III	31 (32%)	33 (32%)
1 n (%)
2 Standardized Mean Difference
3 IC = intervalle de confiance

Table 19.7: différence moyenne standardisée

Pour calculer la différence des proportions pour chaque modalité de grade, il est nécessaire de transformer, en amont, la variable catégorielle grade en trois variables dichotomiques (de type oui/non, une par modalité), ce qui peut se faire facilement avec la fonction fastDummies::dummy_cols() de l’extension fastDummies.

trial |> 
  fastDummies::dummy_cols("grade") |> 
  tbl_summary(
    by = trt,
    include = starts_with("grade_"),
    digits = ~ c(0, 1)
  ) |> 
  add_difference()

Caractéristique	Drug A N = 98 1	Drug B N = 102 1	Difference 2	95% IC 2,3	p-valeur 2
grade_I	35 (35,7%)	33 (32,4%)	3,4%	-11% – 17%	0,7
grade_II	32 (32,7%)	36 (35,3%)	-2,6%	-17% – 11%	0,8
grade_III	31 (31,6%)	33 (32,4%)	-0,72%	-14% – 13%	>0,9
1 n (%)
2 2-sample test for equality of proportions with continuity correction
3 IC = intervalle de confiance

Table 19.8: différence entre proportions avec création de variables dichotomiques

19.2 Une variable continue selon une variable catégorielle

19.2.1 Tableau comparatif avec gtsummary

Dans le chapitre sur la statistique univariée (cf. Section 18.2), nous avons abordé comment afficher les statistiques descriptives d’une variable continue avec gtsummary::tbl_summary(). Pour comparer une variable continue selon plusieurs groupes définis par une variable catégorielle, il suffit d’utiliser le paramètre by :

trial |> 
  tbl_summary(
    include = age,
    by = grade
  )

Caractéristique	I N = 68 1	II N = 68 1	III N = 64 1
Age	47 (37 – 56)	49 (37 – 57)	47 (38 – 58)
Manquant	2	6	3
1 Médiane (Q1 – Q3)

Table 19.9: âge médian et intervalle interquartile selon le grade

La fonction gtsummary::add_overall() permet d’ajouter une colonne total et gtsummary::modify_spanning_header() peut-être utilisé pour ajouter un en-tête de colonne.

trial |> 
  tbl_summary(
    include = age,
    by = grade
  ) |> 
  add_overall(last = TRUE) |> 
  modify_spanning_header(
    all_stat_cols(stat_0 = FALSE) ~ "**Grade**"
  )

Caractéristique	Grade			Overall N = 200 1
	I N = 68 1	II N = 68 1	III N = 64 1
Age	47 (37 – 56)	49 (37 – 57)	47 (38 – 58)	47 (38 – 57)
Manquant	2	6	3	11
1 Médiane (Q1 – Q3)

Table 19.10: âge médian et intervalle interquartile selon le grade

Comme pour un tri à plat, on peut personnaliser les statistiques à afficher avec statistic.

trial |> 
  tbl_summary(
    include = age,
    by = grade,
    statistic = all_continuous() ~ "{mean} ({sd})",
    digits = all_continuous() ~ c(1, 1)
  ) |> 
  add_overall(last = TRUE)

Caractéristique	I N = 68 1	II N = 68 1	III N = 64 1	Overall N = 200 1
Age	46,2 (15,2)	47,5 (13,7)	48,1 (14,1)	47,2 (14,3)
Manquant	2	6	3	11
1 Moyenne (ET)

Table 19.11: âge moyen et écart-type selon le grade

19.2.2 Représentations graphiques

La moyenne ou la médiane sont des indicateurs centraux et ne suffisent pas à rendre compte des différences de distribution d’une variable continue entre plusieurs sous-groupes.

Une représentation usuelle pour comparer deux distributions consiste à avoir recours à des boîtes à moustaches que l’on obtient avec ggplot2::geom_boxplot().

ggplot(trial) +
  aes(x = grade, y = age) +
  geom_boxplot(fill = "lightblue") +
  theme_light()

Figure 19.4: boîtes à moustache

Astuce

Le trait central représente la médiane, le rectangle est délimité par le premier et le troisième quartiles (i.e. le 25^e et le 75^e percentiles). Les traits verticaux vont jusqu’aux extrêmes (minimum et maximum) ou jusqu’à 1,5 fois l’intervalle interquartile. Si des points sont situés à plus d’1,5 fois l’intervalle interquartile au-dessus du 3^e quartile ou en-dessous du 1^er quartile, ils sont considérés comme des valeurs atypiques et représentés par un point. Dans l’exemple précédent, c’est le cas des deux plus petites valeurs observées pour le grade I.

Alternativement, on peut utiliser un graphique en violons qui représentent des courbes de densité dessinées en miroir.

ggplot(trial) +
  aes(x = grade, y = age) +
  geom_violin(fill = "lightblue") +
  theme_light()

Figure 19.5: graphique en violons

Il est toujours possible de représenter les observations individuelles sous la forme d’un nuage de points. Le paramètre alpha permet de rendre les points transparents afin de mieux visualiser les superpositions de points.

ggplot(trial) +
  aes(x = grade, y = age) +
  geom_point(alpha = .25, colour = "blue") +
  theme_light()

Figure 19.6: un nuage de points avec une variable continue et une variable catégorielle

Comme la variable grade est catégorielle, tous les points d’une même modalité sont représentées sur une même ligne. La représentation peut être améliorée en ajoutant un décalage aléatoire sur l’axe horizontal. Cela s’obtient avec ggplot2::position_jitter() en précisant height = 0 pour ne pas ajouter de décalage vertical et width = .2 pour décaler horizontalement les points entre -20% et +20%.

ggplot(trial) +
  aes(x = grade, y = age) +
  geom_point(
    alpha = .25,
    colour = "blue",
    position = position_jitter(height = 0, width = .2)
  ) +
  theme_light()

Figure 19.7: un nuage de points avec une variable continue et une variable catégorielle et avec un décalage horizontal aléatoire

La statistique ggstats::stat_weighted_mean() de ggstats permets de calculer à la volée la moyenne du nuage de points.

ggplot(trial) +
  aes(x = grade, y = age) +
  geom_point(stat = "weighted_mean", colour = "blue") +
  theme_light()

Figure 19.8: âge moyen selon le grade

Cela peut être utile pour effectuer des comparaisons multiples.

ggplot(trial) +
  aes(x = grade, y = age, colour = stage, group = stage) +
  geom_line(stat = "weighted_mean") +
  geom_point(stat = "weighted_mean") +
  facet_grid(cols = vars(trt)) +
  theme_light()

Figure 19.9: âge moyen selon le grade, par traitement et état d’avancement de la maladie

Pyramide des âges

Il est possible de réaliser assez facilement une pyramide des âges en combinant un histogramme avec position_likert_count() fournie par le package ggstats.

Nous allons pour illustrer cela prendre le jeu de données hdv2003 fourni par le package questionr.

data("hdv2003", package = "questionr")
library(ggplot2)
library(ggstats)
ggplot(hdv2003) +
  aes(fill = sexe, y = age) +
  geom_histogram(
    position = "likert_count",
    binwidth = 1,
    color = "black"
  ) +
  scale_x_continuous(label = label_number_abs()) +
  scale_y_continuous(breaks = 1:10 * 10)

19.2.3 Calcul manuel

Le plus simple pour calculer des indicateurs par sous-groupe est d’avoir recours à dplyr::summarise() avec dplyr::group_by().

library(dplyr)
trial |>
  group_by(grade) |> 
  summarise(
    age_moy = mean(age, na.rm = TRUE),
    age_med = median(age, na.rm = TRUE)
  )

# A tibble: 3 × 3
  grade age_moy age_med
  <fct>   <dbl>   <dbl>
1 I        46.2    47  
2 II       47.5    48.5
3 III      48.1    47  

En base R, on peut avoir recours à tapply(). On lui indique d’abord le vecteur sur lequel on souhaite réaliser le calcul, puis un facteur qui indiquera les sous-groupes, puis une fonction qui sera appliquée à chaque sous-groupe et enfin, optionnellement, des arguments additionnels qui seront transmis à cette fonction.

tapply(trial$age, trial$grade, mean, na.rm = TRUE)

       I       II      III 
46.15152 47.53226 48.11475 

19.2.4 Tests de comparaison

Pour comparer des moyennes ou des médianes, le plus facile est encore d’avoir recours à gtsummary et sa fonction gtsummary::add_p().

trial |> 
  tbl_summary(
    include = age,
    by = grade
  ) |> 
  add_p()

Caractéristique	I N = 68 1	II N = 68 1	III N = 64 1	p-valeur 2
Age	47 (37 – 56)	49 (37 – 57)	47 (38 – 58)	0,8
Manquant	2	6	3
1 Médiane (Q1 – Q3)
2 Test de Kruskal-Wallis

Table 19.12: test de comparaison sur la somme des rangs

Par défaut, pour les variables continues, un test de Kruskal-Wallis calculé avec la fonction stats::kruskal.test() est utilisé lorsqu’il y a trois groupes ou plus, et un test de Wilcoxon-Mann-Whitney calculé avec stats::wilcox.test() (test de comparaison des rangs) lorsqu’il n’y a que deux groupes. Au sens strict, il ne s’agit pas de tests de comparaison des médianes mais de tests sur la somme des rangs¹. En pratique, ces tests sont appropriés lorsque l’on présente les médianes et les intervalles inter-quartiles.

¹ Si l’on a besoin spécifiquement d’un test de comparaison des médianes, il existe le test de Brown-Mood disponible dans le package coin avec la fonction coin::median_test(). Attention, il ne faut pas confondre ce test avec le test de dispersion de Mood implémenté dans la fonction stats::mood.test().

Si l’on affiche des moyennes, il serait plus juste d’utiliser un test t de Student (test de comparaison des moyennes) calculé avec stats::t.test(), valable seulement si l’on compare deux moyennes. Pour tester si trois moyennes ou plus sont égales, on aura plutôt recours à stats::oneway.test().

On peut indiquer à gtsummary::add_p() le test à utiliser avec le paramètre test.

trial |> 
  tbl_summary(
    include = age,
    by = grade,
    statistic = all_continuous() ~ "{mean} ({sd})"
  ) |> 
  add_p(
    test = all_continuous() ~ "oneway.test"
  )

Caractéristique	I N = 68 1	II N = 68 1	III N = 64 1	p-valeur 2
Age	46 (15)	48 (14)	48 (14)	0,7
Manquant	2	6	3
1 Moyenne (ET)
2 One-way analysis of means (not assuming equal variances)

Table 19.13: test de comparaison des moyennes

Précision statistique

Classiquement, le test t de Student présuppose l’égalité des variances entre les deux sous-groupes, ce qui permet de former une estimation commune de la variance des deux échantillons (on parle de pooled variance), qui revient à une moyenne pondérée des variances estimées à partir des deux échantillons. Pour tester l’égalité des variances de deux échantillons, on peut utiliser stats::var.test().

Dans le cas où l’on souhaite relaxer cette hypothèse d’égalité des variances, le test de Welch ou la correction de Satterthwaite reposent sur l’idée que l’on utilise les deux estimations de variance séparément, suivie d’une approximation des degrés de liberté pour la somme de ces deux variances.

Par défaut, la fonction stats::t.test() réalise un test de Welch. Pour un test classique de Student, il faut lui préciser var.equal = TRUE.

De manière similaire, stats::oneway.test() ne présuppose pas, par défaut, l’égalité des variances et généralise donc le test de Welch au cas à trois modalités ou plus. Cependant, on peut là encore indiquer var.equal = TRUE, auquel cas une analyse de variance (ANOVA) classique sera réalisée, que l’on peut aussi obtenir avec stats::aov().

Il est possible d’indiquer à gtsummary::add_p() des arguments additionnels à passer à la fonction utilisée pour réaliser le test :

trial |> 
  tbl_summary(
    include = age,
    by = trt,
    statistic = all_continuous() ~ "{mean} ({sd})"
  ) |> 
  add_p(
    test = all_continuous() ~ "t.test",
    test.args = all_continuous() ~ list(var.equal = TRUE)
  )

Caractéristique	Drug A N = 98 1	Drug B N = 102 1	p-valeur 2
Age	47 (15)	47 (14)	0,8
Manquant	7	4
1 Moyenne (ET)
2 Two Sample t-test

19.2.5 Différence de deux moyennes

La fonctions gtsummary::add_difference() permet, pour une variable continue et si la variable catégorielle spécifiée via by n’a que deux modalités, de calculer la différence des deux moyennes, l’intervalle de confiance de cette différence et test si cette différence est significativement différente de 0 avec stats::t.test().

trial |> 
  tbl_summary(
    include = age,
    by = trt,
    statistic = all_continuous() ~ "{mean} ({sd})"
  ) |> 
  add_difference()

Caractéristique	Drug A N = 98 1	Drug B N = 102 1	Difference 2	95% IC 2,3	p-valeur 2
Age	47 (15)	47 (14)	-0,44	-4,6 – 3,7	0,8
Manquant	7	4
1 Moyenne (ET)
2 test de Student
3 IC = intervalle de confiance

Table 19.14: différence de deux moyennes

19.3 Deux variables continues

19.3.1 Représentations graphiques

La comparaison de deux variables continues se fait en premier lieu graphique, en représentant, via un nuage de points, l’ensemble des couples de valeurs. Notez ici l’application d’un niveau de transparence (alpha) afin de faciliter la lecture des points superposés.

ggplot(iris) +
  aes(x = Petal.Length, y = Petal.Width) +
  geom_point(colour = "blue", alpha = .25) +
  theme_light()

Figure 19.10: nuage de points

La géométrie ggplot2::geom_smooth() permets d’ajouter une courbe de tendance au graphique, avec son intervalle de confiance. Par défaut, il s’agit d’une régression polynomiale locale obtenue avec stats::loess().

ggplot(iris) +
  aes(x = Petal.Length, y = Petal.Width) +
  geom_smooth() +
  geom_point(colour = "blue", alpha = .25) +
  theme_light()

Figure 19.11: nuage de points avec une courbe de tendance

Pour afficher plutôt la droite de régression linéaire entre les deux variables, on précisera method = "lm".

ggplot(iris) +
  aes(x = Petal.Length, y = Petal.Width) +
  geom_smooth(method = "lm") +
  geom_point(colour = "blue", alpha = .25) +
  theme_light()

Figure 19.12: nuage de points avec droite de régression linéaire

Astuce pour afficher l’intercept

Supposons que nous souhaitions montrer l’endroit où la droite de régression coupe l’axe des ordonnées (soit le point sur l’axe y pour x = 0).

Nous pouvons étendre la surface du graphique avec ggplot2::expand_limits(). Cependant, cela n’étend pas pour autant la droite de régression.

ggplot(iris) +
  aes(x = Petal.Length, y = Petal.Width) +
  geom_smooth(method = "lm") +
  geom_point(colour = "blue", alpha = .25) +
  theme_light() +
  expand_limits(x = 0, y = -0.5)

`geom_smooth()` using formula = 'y ~ x'

Une solution simple consiste à utiliser l’option fullrange = TRUE dans ggplot2::geom_smooth() pour étendre la droite de régression à l’ensemble du graphique.

ggplot(iris) +
  aes(x = Petal.Length, y = Petal.Width) +
  geom_smooth(method = "lm", fullrange = TRUE) +
  geom_point(colour = "blue", alpha = .25) +
  theme_light() +
  expand_limits(x = 0, y = -0.5)

`geom_smooth()` using formula = 'y ~ x'

On peut contrôler plus finement la partie de droite à afficher avec l’argument xseq (liste des valeurs pour lesquelles on prédit et affiche le lissage). On peut coupler deux appels à ggplot2::geom_smooth() pour afficher l’extension de la droite vers la gauche en pointillés. L’option se = FALSE permet de ne pas calculer d’intervalles de confiance pour ce second appel.

ggplot(iris) +
  aes(x = Petal.Length, y = Petal.Width) +
  geom_smooth(
    method = "lm", 
    xseq = seq(0, 1, by = .1),
    linetype = "dotted",
    se = FALSE
  ) +
  geom_smooth(method = "lm") +
  geom_point(colour = "blue", alpha = .25) +
  theme_light() +
  expand_limits(x = 0, y = -0.5)

`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

La géométrie ggplot2::geom_rug() permet d’afficher une représentation synthétique de la densité de chaque variable sur les deux axes.

ggplot(iris) +
  aes(x = Petal.Length, y = Petal.Width) +
  geom_smooth(method = "lm") +
  geom_point(colour = "blue", alpha = .25) +
  geom_rug() +
  theme_light()

Figure 19.13: nuage de points avec représentation synthétique des densités marginales

19.3.2 Tester la relation entre les deux variables

Si l’on a besoin de calculer le coefficient de corrélation de Pearson entre deux variables, on aura recours à stats::cor().

cor(iris$Petal.Length, iris$Petal.Width)

[1] 0.9628654

Pour aller plus loin, on peut calculer une régression linéaire entre les deux variables avec stats::lm().

m <- lm(Petal.Length ~ Petal.Width, data = iris)
summary(m)

Call:
lm(formula = Petal.Length ~ Petal.Width, data = iris)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.33542 -0.30347 -0.02955  0.25776  1.39453 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  1.08356    0.07297   14.85   <2e-16 ***
Petal.Width  2.22994    0.05140   43.39   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.4782 on 148 degrees of freedom
Multiple R-squared:  0.9271,    Adjusted R-squared:  0.9266 
F-statistic:  1882 on 1 and 148 DF,  p-value: < 2.2e-16

Les résultats montrent une corrélation positive et significative entre les deux variables.

Pour une présentation propre des résultats de la régression linéaire, on utilisera gtsummary::tbl_regression(). La fonction gtsummary::add_glance_source_note() permet d’ajouter différentes statistiques en notes du tableau de résultats.

m |> 
  tbl_regression() |> 
  add_glance_source_note()

Caractéristique	Beta	95% IC 1	p-valeur
Petal.Width	2,2	2,1 – 2,3	<0,001
R² = 0,927; Adjusted R² = 0,927; Sigma = 0,478; Statistique = 1 882; p-valeur = <0,001; df = 1; Log-likelihood = -101; AIC = 208; BIC = 217; Deviance = 33,8; degrés de liberté des résidus = 148; No. Obs. = 150
1 IC = intervalle de confiance

19.4 Matrice de corrélations

Le package GGally et sa fonction GGally::ggpairs() permettent de représenter facilement une matrice de corrélation entre plusieurs variables, tant quantitatives que qualitatives.

library(GGally)
ggpairs(iris)

Figure 19.14: une matrice de corrélation avec ggpairs()

GGally::ggpairs() et sa petite sœur GGally::ggduo() offrent de nombreuses options de personnalisation qui sont détaillées sur le site dédié du package.

ggpairs(trial, mapping = aes(colour = trt))

Figure 19.15: un second example de matrice de corrélation

19.5 webin-R

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

Figure 19.1: un graphique en barres croisant deux variables Figure 19.2: un graphique avec des barres côte à côte Figure 19.3: un graphique en barres cumulées Figure 19.4: boîtes à moustache Figure 19.5: graphique en violons Figure 19.6: un nuage de points avec une variable continue et une variable catégorielle Figure 19.7: un nuage de points avec une variable continue et une variable catégorielle et avec un décalage horizontal aléatoire Figure 19.8: âge moyen selon le grade Figure 19.9: âge moyen selon le grade, par traitement et état d’avancement de la maladie Figure 19.10: nuage de points Figure 19.11: nuage de points avec une courbe de tendance Figure 19.12: nuage de points avec droite de régression linéaire Figure 19.13: nuage de points avec représentation synthétique des densités marginales Figure 19.14: une matrice de corrélation avec ggpairs() Figure 19.15: un second example de matrice de corrélation