Ce chapitre est évoqué dans le webin-R #03 (statistiques descriptives avec gtsummary et esquisse) sur YouTube.

Nous utiliserons dans ce chapitre les données de l’enquête Histoire de vie 2003 fournies avec l’extension questionr.

library(questionr)
data("hdv2003")
d <- hdv2003

Comparaison de moyennes

On peut calculer la moyenne d’âge des deux groupes en utilisant la fonction tapply1 :

tapply(d$age, d$hard.rock, mean)
     Non      Oui 
48.30211 27.57143

L’écart est important. Est-il statistiquement significatif ? Pour cela on peut faire un test t de Student de comparaison de moyennes à l’aide de la fonction t.test :

t.test(age ~ hard.rock, data = d)

    Welch Two Sample t-test

data:  age by hard.rock
t = 9.6404, df = 13.848, p-value = 1.611e-07
alternative hypothesis: true difference in means between group Non and group Oui is not equal to 0
95 percent confidence interval:
 16.11379 25.34758
sample estimates:
mean in group Non mean in group Oui 
         48.30211          27.57143

Le test est extrêmement significatif. L’intervalle de confiance à 95 % de la différence entre les deux moyennes va de 16,1 ans à 25,3 ans.

La valeur affichée pour p est de 1.611e-07. Cette valeur peut paraître étrange pour les non avertis. Cela signifie tout simplement 1,611 multiplié par 10 à la puissance -7, autrement dit 0,0000001611. Cette manière de représenter un nombre est couramment appelée notation scientifique.

Pour plus de détails, voir http://fr.wikipedia.org/wiki/Notation_scientifique.

Il est possible de désactiver la notation scientifique avec la commande :

options(scipen = 999)

Pour rétablir la notation scientifique :

options(scipen = 0)

Nous sommes cependant allés un peu vite en besogne, car nous avons négligé une hypothèse fondamentale du test t : les ensembles de valeur comparés doivent suivre approximativement une loi normale et être de même variance2.

Dans le test de Student, on suppose l’égalité des variances parentes, 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. Dans le cas où l’on souhaite relaxer cette hypothèse, 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. Le même principe s’applique dans le cas de l’analyse de variance à un facteur (cf. oneway.test).

Comment vérifier que l’hypothèse de normalité est acceptable pour ces données ? D’abord avec un petit graphique composés de deux histogrammes :

par(mfrow = c(1, 2))
hist(d$age[d$hard.rock == "Oui"], main = "Hard rock", col = "red")
hist(d$age[d$hard.rock == "Non"], main = "Sans hard rock", col = "red")
Distribution des âges pour appréciation de la normalité

Une alternative consisterait à utiliser des graphiques de type QQ-plot, à l’aide de la fonction qnorm, même si leur utilisation et leur interprétation ne sera pas détaillée ici.

La fonction par permet de modifier de nombreux paramètres graphiques. Ici, l’instruction par(mfrow = c(1, 2)) sert à indiquer que l’on souhaite afficher deux graphiques sur une même fenêtre, plus précisément que la fenêtre doit comporter une ligne et deux colonnes.

Ça a l’air à peu près bon pour les « Sans hard rock », mais un peu plus limite pour les fans de Metallica, dont les effectifs sont d’ailleurs assez faibles. Si on veut en avoir le cœur net on peut utiliser le test de normalité de Shapiro-Wilk avec la fonction shapiro.test :

shapiro.test(d$age[d$hard.rock == "Oui"])

    Shapiro-Wilk normality test

data:  d$age[d$hard.rock == "Oui"]
W = 0.86931, p-value = 0.04104
shapiro.test(d$age[d$hard.rock == "Non"])

    Shapiro-Wilk normality test

data:  d$age[d$hard.rock == "Non"]
W = 0.98141, p-value = 2.079e-15

Visiblement, le test estime que les distributions ne sont pas suffisamment proches de la normalité dans les deux cas.

Et concernant l’égalité des variances ?

tapply(d$age, d$hard.rock, var)
      Non       Oui 
285.62858  62.72527

L’écart n’a pas l’air négligeable. On peut le vérifier avec le test d’égalité des variances fourni par la fonction var.test :

var.test(age ~ hard.rock, data = d)

    F test to compare two variances

data:  age by hard.rock
F = 4.5536, num df = 1985, denom df = 13, p-value =
0.003217
alternative hypothesis: true ratio of variances is not equal to 1
95 percent confidence interval:
 1.751826 8.694405
sample estimates:
ratio of variances 
          4.553644

La différence est très significative. En toute rigueur le test t n’aurait donc pas pu être utilisé. Cela dit, il convient de rappeler que ce test statistique (1) suppose la normalité des distributions et (2) considère comme hypothèse nulle l’égalité des variances (parentes) – ce que l’on souhaiterait vérifier alors qu’on ne peut pas accepter l’hypothèse nulle dans un cadre d’inférence fréquentiste – sans que l’on définisse réellement ce que signifie des variances différentes sur le plan pratique. Est-ce qu’une variation de la variance du simple au double est pertinente au regard du domaine d’étude, ou bien faut-il décider qu’à partir d’un rapport de 4 on peut considérer qu’il y a bien une différence importante entre deux variances ? Sans avoir fixé au préalable cette hypothèse alternative, on ne peut guère conclure à partir de ce test. Une alternative consiste à comparer la forme des distributions à l’aide, par exemple, de diagrammes de type boîtes à moustaches.

Comparaison des rangs ou de la médiane

Damned ! Ces maudits tests statistiques vont-ils nous empêcher de faire connaître au monde entier notre fabuleuse découverte sur l’âge des fans de Sepultura ? Non ! Car voici qu’approche à l’horizon un nouveau test, connu sous le nom de Wilcoxon/Mann-Whitney. Celui-ci a l’avantage d’être non-paramétrique, c’est à dire de ne faire aucune hypothèse sur la distribution des échantillons comparés, à l’exception que celles-ci ont des formes à peu près comparables (essentiellement en termes de variance). Attention, il ne s’agit pas d’un test comparant les différences de médianes (pour cela il existe le test de Mood) mais d’un test reposant sur la somme des rangs des observations, au lieu des valeurs brutes, dans les deux groupes, via la fonction wilcox.test :

wilcox.test(age ~ hard.rock, data = d)

    Wilcoxon rank sum test with continuity correction

data:  age by hard.rock
W = 23980, p-value = 2.856e-06
alternative hypothesis: true location shift is not equal to 0

Ouf ! La différence est hautement significative3. Nous allons donc pouvoir entamer la rédaction de notre article pour la Revue française de sociologie.

Note : le test de Wilcoxon n’est pas adapté pour comparer les rangs lorsque l’on a trois groupes ou plus. On pourra dans ce cas là avoir recours au test de Kruskal-Wallis avec la fonction krukal.test.

kruskal.test(age ~ hard.rock, data = d)

    Kruskal-Wallis rank sum test

data:  age by hard.rock
Kruskal-Wallis chi-squared = 21.913, df = 1, p-value
= 2.852e-06

Note 2 : le test de Mood mentionné plus haut peut être réalisé avec mood.test.

mood.test(age ~ hard.rock, data = d)

    Mood two-sample test of scale

data:  age by hard.rock
Z = -3.1425, p-value = 0.001675
alternative hypothesis: two.sided

Comparaison de proportions

La fonction prop.test, que nous avons déjà rencontrée pour calculer l’intervalle de confiance d’une proportion (voir le chapitre dédié aux intervalles de confiance) permets également d’effectuer un test de comparaison de deux proportions.

Supposons que l’on souhaite comparer la proportion de personnes faisant du sport entre ceux qui lisent des bandes dessinées et les autres :

tab <- xtabs(~ lecture.bd + sport, data = d)
lprop(tab)
          sport
lecture.bd Non   Oui   Total
  Non       64.2  35.8 100.0
  Oui       48.9  51.1 100.0
  Ensemble  63.8  36.1 100.0

Une représentation graphique sous forme de diagramme en barres peut être définie comme suit :

barplot(prop.table(tab, margin = 1) * 100, beside = TRUE, ylim = c(0, 100), xlab = "Sport", legend.text = c("Lecture : non", "Lecture : oui"))
Répartition des individus (en %) selon les variables sport et lecture

Il suffit de transmettre notre tableau croisé (à 2×2 dimensions) à prop.test :

prop.test(tab)

    2-sample test for equality of proportions with
    continuity correction

data:  tab
X-squared = 4, df = 1, p-value = 0.0455
alternative hypothesis: two.sided
95 percent confidence interval:
 -0.002652453  0.308107236
sample estimates:
   prop 1    prop 2 
0.6420891 0.4893617

On pourra également avoir recours à la fonction fisher.test qui renverra notamment l’odds ratio et son intervalle de confiance correspondant :

fisher.test(tab)

    Fisher's Exact Test for Count Data

data:  tab
p-value = 0.0445
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
 1.003372 3.497759
sample estimates:
odds ratio 
  1.871433

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.

On pourra aussi avoir recours à la fonction odds.ratio de l’extension questionr qui réalise le même calcul mais présente le résultat légèrement différemment :

odds.ratio(tab)

Note : pour le calcul du risque relatif, on pourra regarder du côté de la fonction relrisk de l’extension mosaic.

χ² 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 χ² de Pearson4. Celui-ci s’obtient grâce à la fonction chisq.test, appliquée au tableau croisé obtenu avec table ou xtabs5 :

d$qualreg <- as.character(d$qualif)
d$qualreg[d$qualif %in% c("Ouvrier specialise", "Ouvrier qualifie")] <- "Ouvrier"
d$qualreg[d$qualif %in% c(
  "Profession intermediaire",
  "Technicien"
)] <- "Intermediaire"

tab <- table(d$sport, d$qualreg)
tab
     
      Autre Cadre Employe Intermediaire Ouvrier
  Non    38   117     401           127     381
  Oui    20   143     193           119     114
chisq.test(tab)

    Pearson's Chi-squared test

data:  tab
X-squared = 96.798, df = 4, p-value < 2.2e-16

Le test est hautement significatif : on ne peut donc pas considérer qu’il y a indépendance entre les lignes et les colonnes du tableau.

Notons que l’agrégation des niveaux d’une variable catégorielle peut être réalisée d’une manière différente en utilisant les fonctions de gestion des niveaux d’un facteur. Les expressions précédentes sont donc équivalentes à l’approche ci-après, qui ne nécessite pas de convertir d$qualif en chaîne de caractères :

d$qualreg <- d$qualif
levels(d$qualreg)[1:2] <- "Ouvrier"
levels(d$qualreg)[2:3] <- "Intermédiaire"
tab <- table(d$sport, d$qualreg)

On peut affiner l’interprétation du test en déterminant dans quelle cas l’écart à l’indépendance est le plus significatif en utilisant les résidus du test. Ceux-ci sont notamment affichables avec la fonction chisq.residuals de questionr :

chisq.residuals(tab)
     
      Autre Cadre Employe Intermediaire Ouvrier
  Non  0.11 -3.89    0.95         -2.49    3.49
  Oui -0.15  5.23   -1.28          3.35   -4.70

Les cases pour lesquelles l’écart à l’indépendance est significatif ont un résidu dont la valeur est supérieure à 2 ou inférieure à -2 (le fameux nombre 2 issu de la loi normale, au-delà duquel on s’attend à observer au maximum 2,5 % des observations). Ici on constate que la pratique d’un sport est sur-représentée parmi les cadres et, à un niveau un peu moindre, parmi les professions intermédiaires, tandis qu’elle est sous-représentée chez les ouvriers.

Enfin, on peut calculer le coefficient de contingence de Cramer du tableau, qui présente l’avantage de pouvoir être comparé par la suite à celui calculé sur d’autres tableaux croisés. On peut pour cela utiliser la fonction cramer.v de questionr :

cramer.v(tab)
[1] 0.24199

Pour un tableau à 2×2 entrées, comme discuté plus haut, il est également possible de calculer le test exact de Fisher avec la fonction fisher.test. On peut soit lui passer le résultat de table ou xtabs, soit directement les deux variables à croiser.

lprop(table(d$sexe, d$cuisine))
          
           Non   Oui   Total
  Homme     70.0  30.0 100.0
  Femme     44.5  55.5 100.0
  Ensemble  56.0  44.0 100.0
fisher.test(table(d$sexe, d$cuisine))

    Fisher's Exact Test for Count Data

data:  table(d$sexe, d$cuisine)
p-value < 2.2e-16
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
 2.402598 3.513723
sample estimates:
odds ratio 
  2.903253

Le test du χ² de Pearson étant assez robuste quant aux déviations par rapport aux hypothèses d’applications du test (effectifs théoriques tous ≥ 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.

Données pondérées et l’extension survey

Lorsque l’on utilise des données pondérées, on aura recours à l’extension survey6.

Préparons des données d’exemple :

library(survey)
dw <- svydesign(ids = ~1, data = d, weights = ~poids)

Pour comparer deux moyennes à l’aide d’un test t on aura recours à svyttest :

svyttest(age ~ sexe, dw)

    Design-based t-test

data:  age ~ sexe
t = 2.0404, df = 1998, p-value = 0.04144
alternative hypothesis: true difference in mean is not equal to 0
95 percent confidence interval:
 0.08316148 4.19909687
sample estimates:
difference in mean 
          2.141129

Pour le test de Wilcoxon/Mann-Whitney, on pourra avoir recours à svyranktest :

svyranktest(age ~ hard.rock, dw)

    Design-based KruskalWallis test

data:  age ~ hard.rock
t = -11.12, df = 1998, p-value < 2.2e-16
alternative hypothesis: true difference in mean rank score is not equal to 0
sample estimates:
difference in mean rank score 
                   -0.3636859

On ne peut pas utiliser chisq.test directement sur un tableau généré par svytable. Les effectifs étant extrapolés à partir de la pondération, les résultats du test seraient complètement faussés. Si on veut faire un test du χ² sur un tableau croisé pondéré, il faut utiliser svychisq :

rprop(svytable(~ sexe + clso, dw))
          clso
sexe       Oui   Non   Ne sait pas Total
  Homme     51.6  47.0   1.4       100.0
  Femme     43.9  54.8   1.3       100.0
  Ensemble  47.5  51.1   1.4       100.0
svychisq(~ sexe + clso, dw)

    Pearson's X^2: Rao & Scott adjustment

data:  svychisq(~sexe + clso, dw)
F = 3.3331, ndf = 1.9734, ddf = 3944.9024, p-value =
0.03641

L’extension survey ne propose pas de version adaptée du test exact de Fisher. Pour comparer deux proportions, on aura donc recours au test du χ² :

rprop(svytable(~ lecture.bd + sport, dw))
          sport
lecture.bd Non   Oui   Total
  Non       61.0  39.0 100.0
  Oui       46.8  53.2 100.0
  Ensemble  60.7  39.3 100.0
svychisq(~ lecture.bd + sport, dw)

    Pearson's X^2: Rao & Scott adjustment

data:  svychisq(~lecture.bd + sport, dw)
F = 2.6213, ndf = 1, ddf = 1999, p-value = 0.1056

Tests dans les tableaux de gtsummary

Lorsque l’on réalise un tableau croisé avec tbl_summary ou tbl_svysummary de l’extension gtsummary, il est possible d’ajouter des tests de comparaison avec add_p.

library(gtsummary)
theme_gtsummary_language("fr", decimal.mark = ",", big.mark = " ")
Setting theme `language: fr`
d %>%
  tbl_summary(
    include = c("hard.rock", "age", "sport"),
    by = "hard.rock"
  ) %>%
  add_p()
Caractéristique Non, N = 1 9861 Oui, N = 141 p-valeur2
age 48 (35 – 60) 26 (22 – 30) <0,001
sport >0,9
Non 1 268 (64%) 9 (64%)
Oui 718 (36%) 5 (36%)

1 Médiane (EI); n (%)

2 test de Wilcoxon-Mann-Whitney; test du khi-deux d'indépendance

Il est possible de préciser le type de test à utiliser.

d %>%
  tbl_summary(
    include = c("hard.rock", "age"),
    by = "hard.rock"
  ) %>%
  add_p(test = list(all_continuous() ~ "wilcox.test"))
Caractéristique Non, N = 1 9861 Oui, N = 141 p-valeur2
age 48 (35 – 60) 26 (22 – 30) <0,001

1 Médiane (EI)

2 test de Wilcoxon-Mann-Whitney

Cela fonctionne également avec les données pondérées et les plans d’échantillonnage complexe.

dw %>%
  tbl_svysummary(
    include = c("hard.rock", "age", "sport"),
    by = "hard.rock"
  ) %>%
  add_p()
Caractéristique Non, N = 10 981 9631 Oui, N = 89 2641 p-valeur2
age 45 (31 – 60) 23 (19 – 28) <0,001
sport 0,8
Non 6 656 841 (61%) 57 920 (65%)
Oui 4 325 122 (39%) 31 344 (35%)

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

  1. La fonction tapply est présentée plus en détails dans le chapitre Manipulation de données.↩︎

  2. Concernant cette seconde condition, t.test utilise par défaut un test de Welch qui ne suppose pas l’égalité des variances parentes ; il est toutefois possible d’utiliser le test classique de Student en spécifiant l’option var.equal = TRUE.↩︎

  3. Ce test peut également fournir un intervalle de confiance avec l’option conf.int=TRUE.↩︎

  4. On ne donnera pas plus d’indications sur le test du χ² ici. Les personnes désirant une présentation plus détaillée pourront se reporter (attention, séance d’autopromotion !) à la page suivante : http://alea.fr.eu.org/pages/khi2.↩︎

  5. On peut aussi appliquer directement le test en spécifiant les deux variables à croiser via chisq.test(d$qualreg, d$sport).↩︎

  6. Voir le chapitre dédié aux données pondérées.↩︎