shapiro.test(dat$RM1)
Shapiro-Wilk normality test
data: dat$RM1
W = 0.97451, p-value = 0.4777
On ne donne ici que quelques exemples des tests statistiques simples les plus répandus.
En de nombreuses situations (vérification des hypothèses requises en modélisation linéaire, réalisation de tests paramétriques, …), il est nécessaire d’évaluer l’hypothèse de normalité des données, en amont ou en aval des analyses.
Une méthode possible pour évaluer l’hypothèse de normalité pour une variable donnée (pas la seule, et pas nécessairement la meilleure) est de procéder à des tests de normalité. Il existe beaucoup de tests de normalité mais le plus courant est celui de Shapiro-Wilk, implémenté dans une fonction R de base :
shapiro.test(dat$RM1)
Shapiro-Wilk normality test
data: dat$RM1
W = 0.97451, p-value = 0.4777
Conseil : toujours accompagner ces tests de représentations graphiques.
L’intérêt de R en tant que langage de programmation est aussi d’automatiser les tâches répétitives.
On peut obtenir le résultat de tests de normalité pour toutes les variables numériques du jeu de données d’un seul coup en utilisant la fonction apply()
:
## Isoler d'abord les variables numériques :
<- select(dat, where(is.numeric))
num ## Appliquer un test de Shapiro-Wilk sur chaque colonne :
apply(num, MARGIN = 2, FUN = shapiro.test)
$FM1
Shapiro-Wilk normality test
data: newX[, i]
W = 0.95134, p-value = 0.09121
$TM6
Shapiro-Wilk normality test
data: newX[, i]
W = 0.98629, p-value = 0.8803
$RM1
Shapiro-Wilk normality test
data: newX[, i]
W = 0.97451, p-value = 0.4777
En réalité, il est surtout nécessaire de tester la normalité d’une variable au sein de chaque groupe d’individus (par exemple pour les femmes puis pour les hommes).
Pour cela, on peut se servir de la fonction tapply()
:
tapply(dat$FM1, INDEX = dat$Sexe, FUN = shapiro.test)
$F
Shapiro-Wilk normality test
data: X[[i]]
W = 0.87723, p-value = 0.02869
$M
Shapiro-Wilk normality test
data: X[[i]]
W = 0.93782, p-value = 0.1973
L’instruction ci-dessus exécute la fonction shapiro.test()
sur la variable FM1
, à l’intérieur de chaque groupe défini (“indexé”) par le facteur Sexe
.
La fonction à utiliser est la fonction de base var.test()
:
var.test(FM1 ~ Sexe, data = dat)
F test to compare two variances
data: FM1 by Sexe
F = 0.69319, num df = 16, denom df = 20, p-value = 0.4611
alternative hypothesis: true ratio of variances is not equal to 1
95 percent confidence interval:
0.2722086 1.8582994
sample estimates:
ratio of variances
0.6931902
Remarque. Le test réalisé ci-dessus est destiné à trancher entre l’hypothèse nulle \(\mathcal{H}_0 : \sigma_1^2 = \sigma_2^2\), et l’hypothèse alternative \(\mathcal{H}_1 : \sigma_1^2 \neq \sigma_2^2\).
On utilise aussi parfois une version unilatérale du test, destinée à trancher entre l’hypothèse nulle \(\mathcal{H}_0 : \sigma_1^2 \leq \sigma_2^2\), et l’hypothèse alternative \(\mathcal{H}_1 : \sigma_1^2 > \sigma_2^2\). Pour que ce test ait un sens, il faut auparavant s’assurer que les groupes soient numérotés de telle sorte que l’échantillon numéro 1 ait la plus grande variance empirique. Dans le cas contraire, il faut renuméroter les groupes avec la fonction relevel()
, comme illustré ci-dessous.
$Sexe2 <- relevel(dat$Sexe, ref = "M")
datvar.test(FM1 ~ Sexe2, data = dat,
alternative = "greater")
F test to compare two variances
data: FM1 by Sexe2
F = 1.4426, num df = 20, denom df = 16, p-value = 0.2306
alternative hypothesis: true ratio of variances is greater than 1
95 percent confidence interval:
0.6339536 Inf
sample estimates:
ratio of variances
1.442606
La fonction à utiliser est la fonction de base t.test()
:
t.test(FM1 ~ Sexe, var.equal = TRUE, data = dat)
Two Sample t-test
data: FM1 by Sexe
t = -4.3553, df = 36, p-value = 0.0001056
alternative hypothesis: true difference in means between group F and group M is not equal to 0
95 percent confidence interval:
-45.15418 -16.46173
sample estimates:
mean in group F mean in group M
414.2482 445.0562
Cette fonction peut aussi être un moyen détourné d’obtenir les moyennes de chaque groupe (“sample estimates”).
Si l’hypothèse de normalité semble douteuse, il est préférable d’appliquer des tests non paramétriques. Pour l’équivalent non paramétrique du test de Student, la fonction à utiliser est wilcox.test()
:
wilcox.test(FM1 ~ Sexe, data = dat)
Wilcoxon rank sum exact test
data: FM1 by Sexe
W = 60, p-value = 0.0002933
alternative hypothesis: true location shift is not equal to 0
On souhaite comparer la présence d’objets de parure dans l’échantillon de tombes fouillées en fonction du sexe des individus.
La première étape est de construire une table de contingence, présentant les effectifs bruts pour chaque croisement de modalités :
<- table(dat$Sexe, dat$Presence_parure)
contin print(contin)
Non Oui
F 7 12
M 19 4
Comme les effectifs bruts ne sont pas toujours parlants, on peut préférer une table donnant les distributions marginales en pourcentages (ou plutôt en proportion, entre 0 et 1) :
prop.table(contin, margin = 1)
Non Oui
F 0.3684211 0.6315789
M 0.8260870 0.1739130
prop.table(contin, margin = 2)
Non Oui
F 0.2692308 0.7500000
M 0.7307692 0.2500000
Il semble clair que les objets de parure sont associés aux tombes féminines. Pour vérifier cette hypothèse, on procède à un test du \(\chi^2\) :
chisq.test(contin)
Pearson's Chi-squared test with Yates' continuity correction
data: contin
X-squared = 7.4025, df = 1, p-value = 0.006513
En cas de faibles effectifs, ou si les conditions de Cochran ne sont pas satisfaites (elles le sont ici), réaliser un test exact de Fisher sera préférable :
fisher.test(contin)
Fisher's Exact Test for Count Data
data: contin
p-value = 0.003883
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
0.02236343 0.60861101
sample estimates:
odds ratio
0.1302994
La fonction cor()
permet d’obtenir le coefficient de corrélation entre deux variables :
cor(x = dat$RM1, y = dat$FM1, use = "complete.obs")
[1] 0.4782041
Elle renvoie une matrice de corrélation si plus de deux variables sont fournies en argument :
cor(dat[ , c("RM1", "FM1", "TM6")],
use = "complete.obs")
RM1 FM1 TM6
RM1 1.0000000 0.4782041 0.7092951
FM1 0.4782041 1.0000000 0.6153894
TM6 0.7092951 0.6153894 1.0000000
Sur la façon dont cor()
gère les données manquantes pour calculer une matrice de corrélation (argument use
), cf. l’aide de la fonction.
La “significativité” du coefficient de corrélation (test n’ayant pas toujours un grand intérêt pratique) peut être évaluée avec cor.test()
:
cor.test(dat$FM1, dat$RM1)
Pearson's product-moment correlation
data: dat$FM1 and dat$RM1
t = 3.2213, df = 35, p-value = 0.002757
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
0.1824562 0.6945972
sample estimates:
cor
0.4782041