Comparaison de deux moyennes - Tests de Student

Conditions du tests
- X suit une loi normale ou n est grand (>30)
- préciser en option si les variances sont égales var.equal = T (test de Student) ou inégales var.equal = F (test de Welsh)

Exemple : La moyenne d’âge est-elle différente selon qu’on écoute du hard rock ou non ?

Hypothèses
H0 : moyenne d’âge des personnes qui écoutent du hard rock = moyenne d’âge des personnes qui n’écoutent pas de hard rock
H1 : moyenne d’âge des personnes qui écoutent du hard rock \(\ne\) moyenne d’âge des personnes qui n’écoutent pas de hard rock

OU

\({H0 :}\;\overline{Age}_{Hard Rock=Non} = \overline{Age}_{Hard Rock=Oui}\;\) versus \(\;{H1 :}\;\overline{Age}_{Hard Rock=Non} \neq \overline{Age}_{Hard Rock=Oui}\)

Description des données
On calcul les principaux indicateurs de dispersion pour chaque groupe.

d %>% 
  group_by(hard.rock) %>% 
  summarise(
        min=min(age,na.rm = T),
        p25=quantile(age,0.25,na.rm = T),
        median=median(age,na.rm = T),
        p75=quantile(age,0.75,na.rm = T),
        max=max(age,na.rm = T),
        `ecart-type`=sd(age,na.rm = T),
        moyenne=mean(age,na.rm = T),
        n=sum(!is.na(age))
      ) %>% 
  bind_rows(
    d %>% summarise(
        min=min(age,na.rm = T),
        p25=quantile(age,0.25,na.rm = T),
        median=median(age,na.rm = T),
        p75=quantile(age,0.75,na.rm = T),
        max=max(age,na.rm = T),
        `ecart-type`=sd(age,na.rm = T),
        moyenne=mean(age,na.rm = T),
        n=sum(!is.na(age))
      )
    ) %>% 
  mutate(hard.rock=if_else(is.na(hard.rock),"Ensemble",as.character(hard.rock))) %>% 
  kable(digits = 1) %>% 
  kable_styling("striped") 

Explications sur le code
group_by(hard.rock) indique que la suite des opérations doit s’appliquer à chacune des valeurs de hard.rock
bind_rows(...) permet de concaténer des lignes à notre tableau initial. Ici, cela me permet d’ajouter une ligne pour l’ensemble des données
mutate(hard.rock=if_else(is.na(hard.rock),"Ensemble",as.character(hard.rock))) renomme la dernière ligne Ensemble
kable(digits = 1) indique que les variables en format numerique doivent être arrondies un chiffre après la virgule
kable_styling("striped") permet d’avoir un joli tableau en sortie RMarkdown
L’option na.rm = T indique de calculer les statistiques sur le sous-échantillons des données non manquantes. Ici on n’a pas de données manquantes mais c’est souvent le cas. Si vous êtes en présence de données manquantes et que vous n’avaez pas indiqué cette option, R renvoie NA à la place de l’indicateur que vous souhaitiez calculer.

Lecture
L’âge médian des personnes qui écoutent du hard rock est de 25.5 ans alors que l’âge médian des personnes qui n’écoutent pas de hard rock est de 48 ans.
L’âge moyen des personnes qui écoutent du hard rock est de 27.6 ans alors que l’âge moyen des personnes qui n’écoutent pas de hard rock est de 48.3 ans.

Représentation graphique Ci après un exemple de boîte à moustache.

boxplot(age ~ hard.rock, data=d)

Maintenant on effectue le test de Student pour tester si l’âge moyen est différent dans les deux groupes.

t.test(age ~ hard.rock, data=d, var.equal = F, alternative = "two.sided")

## 
##  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

Interprétation La probabilité de rejeter H0 en se trompant est égale à 1.611e-07 soit à 0.0000001611. Ou dit autrement, la p-value est inférieure à 5%, donc on rejette H0, donc on accepte H1.
On rejette l’hypothèse d’égalité des moyenne (H0), ce qui revient à dire que les moyennes sont significativement différentes (H1). Ce qui n’est pas vraiment étonnant puisque les personnes qui écoutent du hard rock sont nettement plus jeunes dans notre échantillon.

Remarques
Ces remarques sont généralisables à l’ensemble des tests présentés dans la suite du document. Pour alléger le document, je me placerai ensuite toujours dans le cas de deux échantillons indépendants et d’un test bilatéral. Mais il est presque toujours possible en modifiant les options d’indiquer s’il s’agit dun test de conformité sur un seul échantillon, ou s’il s’agit de deux échantillons appariés ou encore d’indiquer que l’on souhaite effectuer un test unilatéral.

• Tests unilatéral ou bilatéral
  On aurait pu vouloir tester si les personnes qui n’écoutent pas du hard rock sont plus âgées que les personnes qui en écoutent. C’est ce qu’on appelle un test unilatéral. Les hypothèses deviennent alors les suivantes.

  Test unilatéral
  H0 : la moyenne d’âge des personnes qui n’écoutent pas du hard-rock est <= à la moyenne d’âge des personnes qui écoutent du hard-rock
  H1 : la moyenne d’âge des personnes qui n’écoutent pas du hard-rock est > à la moyenne d’âge des personnes qui écoutent du hard-rock

  OU

  \({H0 :}\;\overline{Age}_{Hard Rock=Non} \leq \overline{Age}_{Hard Rock=Oui}\;\) versus \(\;{H1 :}\;\overline{Age}_{Hard Rock=Non} > \overline{Age}_{Hard Rock=Oui}\)

  La commande à utilser est alors la suivante (on modifie l’option alternative) :
  t.test(age ~ hard.rock, data=d, var.equal = F, alternative = "greater")

• Echantillons appariés
  Les échantillons appariées concernent les cas de mesures répétées sur un même individu. Le calcul de la variance est alors différent et le test est plus puissant. Par exemple, on mesure le taux de glycémie avant et après un effort sportif d’un ensemble d’individu puis on teste si le taux de glycémie moyen a baissé avant et après l’effort. Il faut alors préciser qu’il s’agit d’un échantillon apparié avec l’option paired=TRUE. Si les mesures sont effectuées sur des échantillons distincts, on parle alors d’échantillons indépendants.

• Test de conformité sur un échantillon
  On aurait pu tester si la moyenne d’âge dans notre échantillon était égale à une valeur de référence.
  Exemple : est-ce que la moyenne d’âge des personnes qui n’écoutent pas du hard rock est égale à 30 ?
  H0 : la moyenne d’âge des personnes qui n’écoutent pas du hard rock est égale à 30
  H1 : la moyenne d’âge des personnes qui n’écoutent pas du hard rock est différente de 30
  OU \({H0 :}\;\overline{Age}=30\;\) versus \(\;{H1 :}\;\overline{Age} \ne 30\)

  La commande à utilser est alors la suivante : t.test(d$age[d$hard.rock=="Non"],mu=30,alternative = "two.sided")

Comparaison de k moyennes

iris %>% 
  group_by(Species) %>% 
  summarise(
        min=min(Sepal.Width,na.rm = T),
        p25=quantile(Sepal.Width,0.25,na.rm = T),
        median=median(Sepal.Width,na.rm = T),
        p75=quantile(Sepal.Width,0.75,na.rm = T),
        max=max(Sepal.Width,na.rm = T),
        `ecart-type`=sd(Sepal.Width,na.rm = T),
        moyenne=mean(Sepal.Width,na.rm = T),
        n=sum(!is.na(Sepal.Width))
      ) %>% 
  bind_rows(
    iris %>% summarise(
        min=min(Sepal.Width,na.rm = T),
        p25=quantile(Sepal.Width,0.25,na.rm = T),
        median=median(Sepal.Width,na.rm = T),
        p75=quantile(Sepal.Width,0.75,na.rm = T),
        max=max(Sepal.Width,na.rm = T),
        `ecart-type`=sd(Sepal.Width,na.rm = T),
        moyenne=mean(Sepal.Width,na.rm = T),
        n=sum(!is.na(Sepal.Width))
      )
    ) %>% 
  mutate(Species=if_else(is.na(Species),"Ensemble",as.character(Species))) %>% 
  kable(digits = 1) %>% 
  kable_styling("striped") 

boxplot(Sepal.Width ~ Species, iris)

model = aov(Sepal.Width ~ Species, iris)
summary(model)

##              Df Sum Sq Mean Sq F value Pr(>F)    
## Species       2  11.35   5.672   49.16 <2e-16 ***
## Residuals   147  16.96   0.115                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

qqnorm(model$res)

shapiro.test(model$res)

## 
##  Shapiro-Wilk normality test
## 
## data:  model$res
## W = 0.98948, p-value = 0.323

bartlett.test(model$res, iris$Species)

## 
##  Bartlett test of homogeneity of variances
## 
## data:  model$res and iris$Species
## Bartlett's K-squared = 2.0911, df = 2, p-value = 0.3515

iris %>% 
  group_by(Species) %>% 
  summarise(
        min=min(Sepal.Width,na.rm = T),
        p25=quantile(Sepal.Width,0.25,na.rm = T),
        median=median(Sepal.Width,na.rm = T),
        p75=quantile(Sepal.Width,0.75,na.rm = T),
        max=max(Sepal.Width,na.rm = T),
        `ecart-type`=sd(Sepal.Width,na.rm = T),
        moyenne=mean(Sepal.Width,na.rm = T),
        n=sum(!is.na(Sepal.Width))
      ) %>% 
  bind_rows(
    iris %>% summarise(
        min=min(Sepal.Width,na.rm = T),
        p25=quantile(Sepal.Width,0.25,na.rm = T),
        median=median(Sepal.Width,na.rm = T),
        p75=quantile(Sepal.Width,0.75,na.rm = T),
        max=max(Sepal.Width,na.rm = T),
        `ecart-type`=sd(Sepal.Width,na.rm = T),
        moyenne=mean(Sepal.Width,na.rm = T),
        n=sum(!is.na(Sepal.Width))
      )
    ) %>% 
  mutate(Species=if_else(is.na(Species),"Ensemble",as.character(Species))) %>% 
  kable(digits = 1) %>% 
  kable_styling("striped") 

boxplot(Sepal.Width ~ Species, iris)

model = aov(Sepal.Width ~ Species, iris)
TukeyHSD(model)

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Sepal.Width ~ Species, data = iris)
## 
## $Species
##                        diff         lwr        upr     p adj
## versicolor-setosa    -0.658 -0.81885528 -0.4971447 0.0000000
## virginica-setosa     -0.454 -0.61485528 -0.2931447 0.0000000
## virginica-versicolor  0.204  0.04314472  0.3648553 0.0087802

plot(TukeyHSD(model))

qqnorm(model$res)

shapiro.test(model$res)

## 
##  Shapiro-Wilk normality test
## 
## data:  model$res
## W = 0.98948, p-value = 0.323

bartlett.test(model$res, iris$Species)

## 
##  Bartlett test of homogeneity of variances
## 
## data:  model$res and iris$Species
## Bartlett's K-squared = 2.0911, df = 2, p-value = 0.3515

Comparaison de deux échantillons indépendants : Test de Fisher-Snedecor

On teste si les variances de deux échantillons sont égales.

H0 : \(\sigma^2_1=\sigma^2_2\)
H1 : \(\sigma^2_1\ne\sigma^2_2\)

Conditions Les données issues des échantillons sont gaussiennes

Exemple On va tester si la variance des largeurs de sépal de l’espèce setosa est égale à la variance des largeurs de sépal de l’espèce versicolor.

On décrit les données

iris %>% 
  filter(Species=="setosa" | Species=="versicolor") %>%
  group_by(Species) %>% 
  summarise(
        min=min(Sepal.Width,na.rm = T),
        p25=quantile(Sepal.Width,0.25,na.rm = T),
        median=median(Sepal.Width,na.rm = T),
        p75=quantile(Sepal.Width,0.75,na.rm = T),
        max=max(Sepal.Width,na.rm = T),
        `ecart-type`=sd(Sepal.Width,na.rm = T),
        moyenne=mean(Sepal.Width,na.rm = T),
        n=sum(!is.na(Sepal.Width))
      ) %>% 
  bind_rows(
    iris %>% summarise(
        min=min(Sepal.Width,na.rm = T),
        p25=quantile(Sepal.Width,0.25,na.rm = T),
        median=median(Sepal.Width,na.rm = T),
        p75=quantile(Sepal.Width,0.75,na.rm = T),
        max=max(Sepal.Width,na.rm = T),
        `ecart-type`=sd(Sepal.Width,na.rm = T),
        moyenne=mean(Sepal.Width,na.rm = T),
        n=sum(!is.na(Sepal.Width))
      )
    ) %>% 
  mutate(Species=if_else(is.na(Species),"Ensemble",as.character(Species))) %>% 
  kable(digits = 1) %>% 
  kable_styling("striped") 

Représentation graphique

i <- iris %>% filter(Species=="setosa" | Species=="versicolor")
#i$Species <- droplevels(i$Species)
boxplot(Sepal.Width ~ Species, data=i)

ggplot(data=i, aes(x = Sepal.Length,y = ..density..,color=Species,fill=Species)) + geom_histogram(alpha=0.3)

Test de Fisher

var.test(Sepal.Width ~ Species, data=i)

## 
##  F test to compare two variances
## 
## data:  Sepal.Width by Species
## F = 1.4592, num df = 49, denom df = 49, p-value = 0.1895
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.828080 2.571444
## sample estimates:
## ratio of variances 
##           1.459233

Vérifications des conditions

shapiro.test(iris$Sepal.Width[iris$Species=="versicolor"])

## 
##  Shapiro-Wilk normality test
## 
## data:  iris$Sepal.Width[iris$Species == "versicolor"]
## W = 0.97413, p-value = 0.338

shapiro.test(iris$Sepal.Width[iris$Species=="setosa"])

## 
##  Shapiro-Wilk normality test
## 
## data:  iris$Sepal.Width[iris$Species == "setosa"]
## W = 0.97172, p-value = 0.2715

Interprétation

Le test de Shapiro permet de vérifier la normalité des données issues d’un échantillon.
H0 : les données sont issues d’une population normalement distribuée [homoscédasticité]
H1 : les données ne sont pas issues d’une population normalement distribuée [hétéroscédasticité]

1. La p-valeur associée au test de normalité de Shapiro nous montre que nos deux échantillons suivent bien une loi normale. En effet, les p-valeurs sont égales à 0.338 et 0.2715, elles sont supérieures à 5%, dans les deux cas on accepte H0, ie la normalité des distributions.
   
2. La p-valeur du test de Fisher est égale à 0.1895, elle est supérieure à 5%, on accepte H0. Les variances sont égales dans les deux échantillons. On parle d’homoscédasticité.

Comparaison de k échantillon : Test de Bartlett

Le test de Bartlett permet de vérifier l’égalité des variances des résidus du modèle.
H0 : les variances des résidus du modèles sont identiques dans les 3 groupes d’espèces [homoscédasticité]
H1 : les variances des résidus du modèles sont différentes pour au moins deux espèces [hétéroscédasticité]
OU
H0 [toutes les variances sont égales ] \(\sigma_{1}=\sigma_{2}=...=\sigma_{n}\)
H1 [au moins une variances est différente] \(\exists \; (i,j) \; _{i \ne j} \; tq \; \sigma_{i} \ne \sigma_{j}\)

Conditions
Les données issues des échantillons sont gaussiennes

Exemple
On va tester si la variance des largeurs de sépal est égale pour les 3 espèces de fleurs.

On décrit les données

iris %>% 
  group_by(Species) %>% 
  summarise(
        min=min(Sepal.Width,na.rm = T),
        p25=quantile(Sepal.Width,0.25,na.rm = T),
        median=median(Sepal.Width,na.rm = T),
        p75=quantile(Sepal.Width,0.75,na.rm = T),
        max=max(Sepal.Width,na.rm = T),
        `ecart-type`=sd(Sepal.Width,na.rm = T),
        moyenne=mean(Sepal.Width,na.rm = T),
        n=sum(!is.na(Sepal.Width))
      ) %>% 
  bind_rows(
    iris %>% summarise(
        min=min(Sepal.Width,na.rm = T),
        p25=quantile(Sepal.Width,0.25,na.rm = T),
        median=median(Sepal.Width,na.rm = T),
        p75=quantile(Sepal.Width,0.75,na.rm = T),
        max=max(Sepal.Width,na.rm = T),
        `ecart-type`=sd(Sepal.Width,na.rm = T),
        moyenne=mean(Sepal.Width,na.rm = T),
        n=sum(!is.na(Sepal.Width))
      )
    ) %>% 
  mutate(Species=if_else(is.na(Species),"Ensemble",as.character(Species))) %>% 
  kable(digits = 1) %>% 
  kable_styling("striped") 

Représentation graphique

boxplot(Sepal.Width ~ Species, data=iris)

ggplot(data=iris, aes(x = Sepal.Length,y = ..density..,color=Species,fill=Species)) + geom_histogram(alpha=0.3)

Test de Bartlett

bartlett.test(Sepal.Width ~ Species, iris)

## 
##  Bartlett test of homogeneity of variances
## 
## data:  Sepal.Width by Species
## Bartlett's K-squared = 2.0911, df = 2, p-value = 0.3515

Vérifications des conditions

shapiro.test(iris$Sepal.Width[iris$Species=="versicolor"])

## 
##  Shapiro-Wilk normality test
## 
## data:  iris$Sepal.Width[iris$Species == "versicolor"]
## W = 0.97413, p-value = 0.338

shapiro.test(iris$Sepal.Width[iris$Species=="setosa"])

## 
##  Shapiro-Wilk normality test
## 
## data:  iris$Sepal.Width[iris$Species == "setosa"]
## W = 0.97172, p-value = 0.2715

shapiro.test(iris$Sepal.Width[iris$Species=="virginica"])

## 
##  Shapiro-Wilk normality test
## 
## data:  iris$Sepal.Width[iris$Species == "virginica"]
## W = 0.96739, p-value = 0.1809

Interprétation

Le test de Shapiro permet de vérifier la normalité des données issues d’un échantillon.
H0 : les données sont issues d’une population normalement distribuée [homoscédasticité]
H1 : les données ne sont pas issues d’une population normalement distribuée [hétéroscédasticité]

1. La p-valeur associée au test de normalité de Shapiro nous montre que nos trois échantillons suivent bien une loi normale. En effet, les p-valeurs sont égales à 0.338, 0.2715 et 0.1809, elles sont supérieures à 5%, dans les trois cas on accepte H0, ie la normalité des distributions.
   
2. La p-valeur du test de Bartlett sur les résidus du modèle est égale à 0.3515. La p-valeur est supérieure à 5%, on accepte donc H0. Les résidus du modèle sont homoscédastiques.

Comparaison de 2 médianes - Wilcoxon

Exemple : L’âge médian est-il différente selon qu’on écoute du hard rock ou non ?

Hypothèses
H0 : age médian des personnes qui écoutent du hard rock = age médian des personnes qui n’écoutent pas de hard rock
H1 : age médian des personnes qui écoutent du hard rock \(\ne\) age médian des personnes qui n’écoutent pas de hard rock

OU

\({H0 :}\;MedianeAge_{Hard Rock=Non} = MedianeAge_{Hard Rock=Oui}\;\) versus \(\;{H1 :}\;MedianeAge_{Hard Rock=Non} \neq MedianeAge_{Hard Rock=Oui}\)

Description des données
On calcul les principaux indicateurs de dispersion pour chaque groupe.

d %>% 
  group_by(hard.rock) %>% 
  summarise(
        min=min(age,na.rm = T),
        p25=quantile(age,0.25,na.rm = T),
        median=median(age,na.rm = T),
        p75=quantile(age,0.75,na.rm = T),
        max=max(age,na.rm = T),
        `ecart-type`=sd(age,na.rm = T),
        moyenne=mean(age,na.rm = T),
        n=sum(!is.na(age))
      ) %>% 
  bind_rows(
    d %>% summarise(
        min=min(age,na.rm = T),
        p25=quantile(age,0.25,na.rm = T),
        median=median(age,na.rm = T),
        p75=quantile(age,0.75,na.rm = T),
        max=max(age,na.rm = T),
        `ecart-type`=sd(age,na.rm = T),
        moyenne=mean(age,na.rm = T),
        n=sum(!is.na(age))
      )
    ) %>% 
  mutate(hard.rock=if_else(is.na(hard.rock),"Ensemble",as.character(hard.rock))) %>% 
  kable(digits = 1) %>% 
  kable_styling("striped") 

Représentation graphique Ci après un exemple de boîte à moustache.

boxplot(age ~ hard.rock, data=d)

Maintenant on effectue le test de Wilcoxon pour tester si l’âge médian est différent dans les deux groupes.

wilcox.test(age ~ hard.rock, data=d, alternative = "two.sided")

## 
##  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

Interprétation

La probabilité de rejeter H0 en se trompant est égale à 2.856e-06 soit à 0.000002856. Ou dit autrement, la p-value est inférieure à 5%, donc on rejette H0, donc on accepte H1.
On rejette l’hypothèse d’égalité des médiane (H0), ce qui revient à dire que les médianes sont significativement différentes (H1).

Comparaison de k médianes Kruskall Wallis

iris %>% 
  group_by(Species) %>% 
  summarise(
        min=min(Sepal.Width,na.rm = T),
        p25=quantile(Sepal.Width,0.25,na.rm = T),
        median=median(Sepal.Width,na.rm = T),
        p75=quantile(Sepal.Width,0.75,na.rm = T),
        max=max(Sepal.Width,na.rm = T),
        `ecart-type`=sd(Sepal.Width,na.rm = T),
        moyenne=mean(Sepal.Width,na.rm = T),
        n=sum(!is.na(Sepal.Width))
      ) %>% 
  bind_rows(
    iris %>% summarise(
        min=min(Sepal.Width,na.rm = T),
        p25=quantile(Sepal.Width,0.25,na.rm = T),
        median=median(Sepal.Width,na.rm = T),
        p75=quantile(Sepal.Width,0.75,na.rm = T),
        max=max(Sepal.Width,na.rm = T),
        `ecart-type`=sd(Sepal.Width,na.rm = T),
        moyenne=mean(Sepal.Width,na.rm = T),
        n=sum(!is.na(Sepal.Width))
      )
    ) %>% 
  mutate(Species=if_else(is.na(Species),"Ensemble",as.character(Species))) %>% 
  kable(digits = 1) %>% 
  kable_styling("striped") 

boxplot(Sepal.Width ~ Species, iris)

kruskal.test(Sepal.Width ~ Species, data=iris)

## 
##  Kruskal-Wallis rank sum test
## 
## data:  Sepal.Width by Species
## Kruskal-Wallis chi-squared = 63.571, df = 2, p-value = 1.569e-14

iris %>% 
  group_by(Species) %>% 
  summarise(
        min=min(Sepal.Width,na.rm = T),
        p25=quantile(Sepal.Width,0.25,na.rm = T),
        median=median(Sepal.Width,na.rm = T),
        p75=quantile(Sepal.Width,0.75,na.rm = T),
        max=max(Sepal.Width,na.rm = T),
        `ecart-type`=sd(Sepal.Width,na.rm = T),
        moyenne=mean(Sepal.Width,na.rm = T),
        n=sum(!is.na(Sepal.Width))
      ) %>% 
  bind_rows(
    iris %>% summarise(
        min=min(Sepal.Width,na.rm = T),
        p25=quantile(Sepal.Width,0.25,na.rm = T),
        median=median(Sepal.Width,na.rm = T),
        p75=quantile(Sepal.Width,0.75,na.rm = T),
        max=max(Sepal.Width,na.rm = T),
        `ecart-type`=sd(Sepal.Width,na.rm = T),
        moyenne=mean(Sepal.Width,na.rm = T),
        n=sum(!is.na(Sepal.Width))
      )
    ) %>% 
  mutate(Species=if_else(is.na(Species),"Ensemble",as.character(Species))) %>% 
  kable(digits = 1) %>% 
  kable_styling("striped") 

boxplot(Sepal.Width ~ Species, iris)

pairwise.wilcox.test(iris$Sepal.Width, iris$Species, p.adjust = "BH")

## 
##  Pairwise comparisons using Wilcoxon rank sum test with continuity correction 
## 
## data:  iris$Sepal.Width and iris$Species 
## 
##            setosa  versicolor
## versicolor 6.4e-13 -         
## virginica  1.1e-08 0.0046    
## 
## P value adjustment method: BH

Test de corrélation : Pearson - test paramétrique

On teste la corrélation linéaire entre deux variables.
H0 : pas de corrélation entre les variables, \(r=0\)
H1 : corrélation linéaire entre les variables, \(r\ne0\)

Conditions
normalité des variables et lien linéaire

Exemple
On va tester la corrélation linéaire entre la longueur et la largeur des pétales.

Description des données

library(skimr)

iris %>% 
  summarise(
        nomVar="Petal.Width",
        min=min(Petal.Width,na.rm = T),
        p25=quantile(Petal.Width,0.25,na.rm = T),
        median=median(Petal.Width,na.rm = T),
        p75=quantile(Petal.Width,0.75,na.rm = T),
        max=max(Petal.Width,na.rm = T),
        `ecart-type`=sd(Petal.Width,na.rm = T),
        moyenne=mean(Petal.Width,na.rm = T),
        n=sum(!is.na(Petal.Width))) %>%
kable(digits = 2) %>% 
kable_styling("striped") 

iris %>% 
summarise(
      nomVar="Petal.Length",
      min=min(Petal.Length,na.rm = T),
      p25=quantile(Petal.Length,0.25,na.rm = T),
      median=median(Petal.Length,na.rm = T),
      p75=quantile(Petal.Length,0.75,na.rm = T),
      max=max(Petal.Length,na.rm = T),
      `ecart-type`=sd(Petal.Length,na.rm = T),
      moyenne=mean(Petal.Length,na.rm = T),
      n=sum(!is.na(Petal.Length))) %>%
kable(digits = 2) %>% 
kable_styling("striped") 

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

[1] 0.9628654

Représentation graphique

plot(iris$Petal.Length, iris$Petal.Width, ylab = "Longueur des pétales", xlab = "Largeur des pétales")
abline(lm(Petal.Width ~ Petal.Length, data = iris), col = "red")

Test de Pearson

cor.test(iris$Petal.Width, iris$Petal.Length, method = "pearson")

## 
##  Pearson's product-moment correlation
## 
## data:  iris$Petal.Width and iris$Petal.Length
## t = 43.387, df = 148, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.9490525 0.9729853
## sample estimates:
##       cor 
## 0.9628654

Interprétation La p-valeur est égale à < 2.2e-16. La probabilité est tellement faible que le logiciel vous indique qu’elle est inférieure à 0.0000000000000002, ce qui s’écrit bien \(2*10^{-16}\) en écriture scientifique. La p-valeur est inférieure à 5%, on rejette H0 et on accepte H1. Il y a bien une corrélation linéaire entre la longueur et la lergeur des pétales.

Test de corrélation : Spearman - test non paramétrique

On teste la corrélation entre les rangs des valeurs.

H0 : pas de corrélation entre les variables, \(\rho=0\)
H1 : corrélation linéaire entre les variables, \(\rho\ne0\)

Exemple On va tester la corrélation entre la longueur et la largeur des pétales.

Description des données

library(skimr)

iris %>% 
  summarise(
        nomVar="Petal.Width",
        min=min(Petal.Width,na.rm = T),
        p25=quantile(Petal.Width,0.25,na.rm = T),
        median=median(Petal.Width,na.rm = T),
        p75=quantile(Petal.Width,0.75,na.rm = T),
        max=max(Petal.Width,na.rm = T),
        `ecart-type`=sd(Petal.Width,na.rm = T),
        moyenne=mean(Petal.Width,na.rm = T),
        n=sum(!is.na(Petal.Width))) %>%
kable(digits = 2) %>% 
kable_styling("striped") 

iris %>% 
summarise(
      nomVar="Petal.Length",
      min=min(Petal.Length,na.rm = T),
      p25=quantile(Petal.Length,0.25,na.rm = T),
      median=median(Petal.Length,na.rm = T),
      p75=quantile(Petal.Length,0.75,na.rm = T),
      max=max(Petal.Length,na.rm = T),
      `ecart-type`=sd(Petal.Length,na.rm = T),
      moyenne=mean(Petal.Length,na.rm = T),
      n=sum(!is.na(Petal.Length))) %>%
kable(digits = 2) %>% 
kable_styling("striped") 

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

[1] 0.9628654

Représentation graphique

plot(iris$Petal.Length, iris$Petal.Width, ylab = "Longueur des pétales", xlab = "Largeur des pétales")
abline(lm(Petal.Width ~ Petal.Length, data = iris), col = "red")

Test de Spearman

cor.test(iris$Petal.Width, iris$Petal.Length, method = "spearman")

## 
##  Spearman's rank correlation rho
## 
## data:  iris$Petal.Width and iris$Petal.Length
## S = 35061, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.9376668

Interprétation La p-valeur est égale à < 2.2e-16. La probabilité est tellement faible que le logiciel vous indique qu’elle est inférieure à 0.0000000000000002, ce qui s’écrit bien \(2*10^{-16}\) en écriture scientifique. La p-valeur est inférieure à 5%, on rejette H0 et on accepte H1. Il y a bien une corrélation entre la longueur et la largeur des pétales.

Test du chi-deux

On teste l’indépendance de deux variables.
H0 : les deux variables sont indépendantes
H1 : les deux variables ne sont pas indépendantes

Conditions
Effectifs théoriques >5

Exemple
On veut savoir si d’être un homme ou une femme est associé au fait d’aller au cinéma.

Tableau croisé des effectifs

library(janitor)
d %>% 
  tabyl(sexe,cinema) %>%
  adorn_totals(c("row", "col")) %>%
  adorn_title("combined")

##  sexe/cinema  Non Oui Total
##        Homme  542 357   899
##        Femme  632 469  1101
##        Total 1174 826  2000

Lecture : Parmi les 1101 femmes, 469 vont au cinéma.

Tableau des pourcentages ligne

library(janitor)
d %>% 
  tabyl(sexe,cinema) %>%
  adorn_percentages("row") %>%
  adorn_totals(c("col")) %>%
  adorn_pct_formatting(rounding = "half up", digits = 1) %>%
  adorn_title("combined")

##  sexe/cinema   Non   Oui  Total
##        Homme 60.3% 39.7% 100.0%
##        Femme 57.4% 42.6% 100.0%

Lecture : Parmi les hommes, 39.7% vont au cinéma.

Représentation graphique

ggplot(d) +
 aes(x = sexe, fill = cinema) +
 geom_bar(position = "fill") +
 scale_fill_hue() +
 labs(y = "Pourcentage", title = "Sortie cinéma selon le sexe") +
 theme_bw()

Tests du chi-deux

chisq.test(table(d$sexe,d$cinema))

## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  table(d$sexe, d$cinema)
## X-squared = 1.5843, df = 1, p-value = 0.2081

Interprétation

La p-valeur est égale à 0.2081, la p-valeur est supérieure à 5%. On accepte HO. Les deux variables sont indépendantes. Les hommes ne vont pas plus au cinéma que les femmes et inversement.

Test exact de Fisher

On teste l’indépendance de deux variables.
H0 : les deux variables sont indépendantes
H1 : les deux variables ne sont pas indépendantes

Le test de Fischer est utilisé sur les échantillons de petites tailles quand les conditions de réalisation du Chi-Deux ne sont par réunies.

Exemple On veut savoir si d’être un homme ou une femme est associé au fait d’aller au cinéma.

Tableau croisé des effectifs

library(janitor)
d %>% 
  tabyl(sexe,cinema) %>%
  adorn_totals(c("row", "col")) %>%
  adorn_title("combined")

##  sexe/cinema  Non Oui Total
##        Homme  542 357   899
##        Femme  632 469  1101
##        Total 1174 826  2000

Lecture : Parmi les 1101 femmes, 469 vont au cinéma.

Tableau des pourcentages ligne

library(janitor)
d %>% 
  tabyl(sexe,cinema) %>%
  adorn_percentages("row") %>%
  adorn_totals(c("col")) %>%
  adorn_pct_formatting(rounding = "half up", digits = 1) %>%
  adorn_title("combined")

##  sexe/cinema   Non   Oui  Total
##        Homme 60.3% 39.7% 100.0%
##        Femme 57.4% 42.6% 100.0%

Lecture : Parmi les hommes, 39.7% vont au cinéma.

Représentation graphique

ggplot(d) +
 aes(x = sexe, fill = cinema) +
 geom_bar(position = "fill") +
 scale_fill_hue() +
 labs(y = "Pourcentage", title = "Sortie cinéma selon le sexe") +
 theme_bw()

Tests exact de Fischer

fisher.test(table(d$sexe,d$cinema))

## 
##  Fisher's Exact Test for Count Data
## 
## data:  table(d$sexe, d$cinema)
## p-value = 0.2013
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
##  0.938038 1.353438
## sample estimates:
## odds ratio 
##   1.126576

Interprétation

La p-valeur est égale à 0.2013, la p-valeur est supérieure à 5%. On accepte H0. Les deux variables sont indépendantes. Les hommes ne vont pas plus au cinéma que les femmes et inversement.

Test de normalité de Shapiro-Wilk

Le test de Shapiro permet de tester si les données issues d’une population sont distribuées normalement

H0 : les données sont issues d’une population normalement distribuée
H1 : les données ne sont pas issues d’une population normalement distribuée

Exemple
On souhaite tester si la ditribution de la largeur des sépales suit une loi normale.

Représentationn graphique

hist(iris$Sepal.Width, prob=T,ylim=c(0,1.2))
curve(dnorm(x,mean(iris$Sepal.Width),sd(iris$Sepal.Width)),add=T,col="red")

Test de Shapiro-Wilk

shapiro.test(iris$Sepal.Width)

## 
##  Shapiro-Wilk normality test
## 
## data:  iris$Sepal.Width
## W = 0.98492, p-value = 0.1012

Interprétation

La p-valeur est égale à 0.1012, elle est supérieure à 5%, on accepte H0. Les données sont normalement distribuées.

Test de Kolmogorov - Smirnov

Il s’agit d’un test non paramétrique qui compare 2 distributions au moyen de la fonction de répartition.

H0 : les deux fonctions de répartitions sont égales
H1 : les deux fonctions de répartition sont différentes

Représentation graphique

On peut représenter la densité de chacune des variables.

ggplot() + geom_histogram(data=iris,aes(x = Sepal.Length,y = ..density..),color = "black", fill = "blue",alpha=0.3) +geom_histogram(data=iris,aes(x = Sepal.Width,y = ..density..),color = "black", fill = "green", alpha=0.3)

On peut aussi représenter les fonction de répartition empirique pour chacune des deux variables.

ggplot() + stat_ecdf(data=iris,aes(x = Petal.Length), geom = "step",colour = "red") + stat_ecdf(data=iris,aes(x = Sepal.Length), geom = "step",colour = "blue")

Test

ks.test(iris$Petal.Length, iris$Sepal.Length)

## 
##  Asymptotic two-sample Kolmogorov-Smirnov test
## 
## data:  iris$Petal.Length and iris$Sepal.Length
## D = 0.56, p-value < 2.2e-16
## alternative hypothesis: two-sided

Interprétation

La p-valeur est égale à < 2.2e-16. La probabilité est tellement faible que le logiciel vous indique qu’elle est inférieure à 0.00000000000000022, ce qui s’écrit bien \(2.2*10^{-16}\) en écriture scientifique. La p-valeur est inférieure à 5%, on rejette H0 et on accepte H1. La distribution des valeurs des longueurs de pétales est différente de la distribution des longueurs des sépales.

Note
On peut aussi comparer la distribution d’une variable à une distribution théorique. Est-ce que la largeur des sépales suit une loi normale ? On utilise la commande suivante :ks.test(iris$Sepal.Width,"pnorm",mean(iris$Sepal.Width),sd(iris$Sepal.Width))
Il s’agit du même exemple que celui du test de Shapiro mais cette fois ci on utilise un test non paramétrique (plus concervateur).