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 comparaison de moyennes à l’aide de la fonction t.test :

t.test(d$age ~ d$hard.rock)

    Welch Two Sample t-test

data:  d$age by d$hard.rock
t = 9.6404, df = 13.848, p-value = 1.611e-07
alternative hypothesis: true difference in means 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 14,5 ans à 21,8 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.

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. Comment le vérifier ?

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é

La fonction par permet de modifier de nombreux paramètres graphiques. 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 coeur 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(d$age ~ d$hard.rock)

    F test to compare two variances

data:  d$age by d$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é.

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. Par contre il ne compare pas des différences de moyennes mais des différences de médianes :

wilcox.test(d$age ~ d$hard.rock)

    Wilcoxon rank sum test with continuity
    correction

data:  d$age by d$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.

Comparaison de proportions

La fonction prop.test, que nous avons déjà rencontrer 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, 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

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(table(d$lecture.bd, d$sport))

    Fisher's Exact Test for Count Data

data:  table(d$lecture.bd, d$sport)
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 

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)
                  OR  2.5 % 97.5 %      p  
Fisher's test 1.8714 1.0034 3.4978 0.0445 *
---
Signif. codes:  
0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

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ée, on peut tester l’existence d’un lien entre les modalités de deux variables, avec le très classique test du χ²4. 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 pas considérer qu’il y a indépendance entre les lignes et les colonnes du tableau.

On peut affiner l’interprétation du test en déterminant dans quelle case 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. 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 sousreprésentée chez les ouvriers.

Enfin, on peut calculer le coefficient de contingence de Cramer du tableau, qui peut nous permettre de le comparer par la suite à 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, il est 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 

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
sample estimates:
difference in mean 
          2.141129 

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.9000,
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

  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 propose une option nommée var.equal qui permet d’utiliser une approximation dans le cas où les variances ne sont pas égales.

  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.