07-tests.Rmd

---
output: html_document
editor_options: 
  chunk_output_type: console
---
# Tests

Les tests sortent du champ de la statistique descriptive "pure" pour entrer dans celui de la statistique inférentielle. L'inférence statistique consiste à induire les caractéristiques d'une population à partir de celles d'un échantillon de cette population.

## Sur une variable quantitative

Plusieurs questions peuvent se poser à propos de la distribution d'une variable quantitative. La forme de la distribution (normale ou pas), ainsi que la tendance centrale peuvent être étudiées.

### La distribution est-elle normale ?

Plusieurs approches sont envisageables pour se faire une idée de la réponse à cette question. 
La première est de visualiser graphiquement la distribution, et son "écart" à une loi normale.

Le diagramme quantile-quantile sert à comparer les quantiles de la variable à ceux d'une distribution attendue, par exemple la loi normale. La fonction `stat_qq()` permet d'afficher ce graphique.

```{r, warning = FALSE, fig.height = 4.5, fig.width = 8}
ggplot (data = dat, aes (sample = log_dens)) +
  stat_qq ()
```

Plus la courbe est rectiligne, plus cela signifie que les quantiles de la variable sont similaires à ceux de la distribution attendue.

On peut aussi apprécier la normalité en repartant de l'histogramme :

```{r, warning = FALSE, fig.height = 4.5, fig.width = 8}
vec_logd <- pull (dat, log_dens)
n <- sum (!is.na (vec_logd)) # nombre d'éléments non-NA dans vec_logd
moy <- mean (vec_logd, na.rm = T)
ecart_t <- sd (vec_logd, na.rm = T)
n_binwidth <-  0.2 

ggplot (dat, aes (x = log_dens)) +
  geom_histogram (binwidth = n_binwidth, 
                 colour = "white", fill = "cornflowerblue", size = 0.1) +
  stat_function (fun = function(x) dnorm (x, mean = moy, sd = ecart_t) * n * n_binwidth,
                color = "darkred", size = 1) +
  xlab (label = "Log (densité)") +
  ylab (label = "Nombre de communes")

```

`stat_function()` trace la fonction indiquée. Ici, on trace une loi normale (`dnorm()`) de même moyenne et écart-type que la variable dont on veut tester la normalité.


Après une visualisation, il est possible de "tester" l'affirmation. On se place dans l'hypothèse la plus simple (appelée **hypothése nulle**, ici : c'est que la distribution est bien normale), et on regarde si cela semble cohérent avec nos données.

Tests de normalité : Shapiro-Wilk, préconisé par [Razali et al. 2011 Journal of Statistical Modeling and Analytics 2 (1): 21–33](https://www.google.com/url?sa=t&rct=j&q=&esrc=s&source=web&cd=1&ved=2ahUKEwjQ5PjRitvcAhWNxoUKHTI7Ca4QFjAAegQIABAC&url=http%3A%2F%2Fwww.de.ufpb.br%2F~ulisses%2Fdisciplinas%2Fnormality_tests_comparison.pdf&usg=AOvVaw3O6RzpOVVa3-kzECx4bBDJ)

```{r, warning = FALSE}
dat$log_dens[1:4000] %>% # maxi 5000 données pour ce test  
  shapiro.test () 

```
p < 5% donc l'hypothèse de normalité est rejetée.

### La moyenne diffère-t-elle de la valeur attendue ?

De la même manière, si on s'intéresse à la moyenne d'une distribution, on "teste" le fait qu'elle différe ou non d'une valeur fixe.
Le test de Student (fonction `t.test()`) pour échantillon unique permet de comparer la moyenne observée à une valeur de référence. Il est utilisable  si et seulement si la distribution des valeurs dans l'échantillon est normale. Si ce n'est pas le cas, il faut utiliser le test de Wilcoxon (sur les rangs) avec la fonction `wilcox.test()`.

Exemple : La densité (log-transformée) en Corse diffère-t-elle de la moyenne nationale ?

```{r, warning = FALSE}
dens_corse <- dat %>% filter (REG == "94") %>% pull(log_dens)
dens_nat <- dat %>% pull(log_dens)
moy_dens_nat <- dens_nat %>% mean(na.rm = T) # la valeur de référence

shapiro.test (dens_corse) # non normal -> Wilcoxon
wilcox.test (dens_corse, mu = moy_dens_nat)
t.test (dens_corse, mu = moy_dens_nat) # à titre indicatif
```

La moyenne observée en Corse est de `r round (mean (dens_corse, na.rm = T), 2)`, contre `r round (mean (dens_nat, na.rm = T), 2)` au niveau national. Cette différence est significative (p-value << 5%).





## Sur une variable qualitative

On travaille cette fois sur une variable qualitative. On peut s'interroger sur la répartition des observations.

### La distribution des observations diffère-t-elle de celle attendue ?

Par exemple sur dans un échantillon de personnes, on veut savoir si le sex-ratio est significativement différent de celui de la population française dans son ensemble, soit 51.4%, ou bien si au contraire le pourcentage observé est dans l'intervalle d'incertitude compte tenu de la taille de l'échantillon.

```{r}
femmes <- 601
hommes <- 514

observed_n <- c(femmes, hommes)        # Observation
expected_p <- c(0.514, 1-0.514)        # Proportions attendues
```

Ici l'échantillon est contitué de `r femmes` femmes et `r hommes` hommes donc un total de `r hommes+femmes` personnes.

```{r, fig.width = 3}
expected_n <- round (expected_p * (femmes + hommes))
tab <- cbind (Sexe = c('Femmes', 'Hommes'), `Observé` = observed_n, Attendu = expected_n) %>%
  as.data.frame () %>%
  datatable ()
tab
```

En appliquant le pourcentage théorique, on s'attendrait à ce que ces `r hommes+femmes` personnes soient réparties en `r expected_n[1]` femmes et `r expected_n[2]` hommes. 

Question : Est-ce que les proportions observées peuvent être dues au hasard, ou bien notre échantillon est-il différent de la population ?
On va tester au moyen du test du $\chi^2$ pour échantillon unique. Ce test s'appelle grâce à la fonction `chisq.test()`.
L'hypothèse nulle est que l'échantillon provient de la population générale.

```{r chi2_one_sample}
test <- chisq.test(x = observed_n, p = expected_p)
test
```

Les proportions observées sont dans l'intervalle d'incertitude pour un échantillon de `r hommes+femmes` personnes (p = `r round (test$p.value, 3)`> 5%). L'hypothèse nulle ne peut donc pas être rejetée au risque de 5%.




## Sur 2 variables qualitatives 

On étudie désormais 2 variables qualitatives.

### Les deux variables qualitatives sont-elles indépendantes l'une de l'autre ?

Pour mesurer le lien entre 2 variables qualitatives, on peut se servir du $\chi^2$ :
- Le $\chi^2$ de Pearson constitue la base de toute mesure de liaison entre des variables qualitatives
$\chi^2 = \sum_{i,j} \frac{(N_{ij}-N_{ij}^{*})^2}{N_{ij}^{*}} \hookrightarrow \chi_{(k-1)(r-1)}^2$ où $N_{ij}^{*}=\frac{N_{i.} \cdot N_{.j}}{N}$
- Il mesure l'écart par rapport à la situation d'indépendance : plus le $\chi^2$ est élevé, plus les variables sont liées (elles ne sont pas indépendantes)
- La comparaison de la valeur avec le fractile d'ordre $1-\alpha$ d'une loi $\chi_{(k-1)(r-1)}^2$ permet de connaître la significativité du lien
- On le réalise sur le tableau de contingence des <strong>effectifs</strong>
- *Il faut un effectif > 5  dans chaque case*

[Cette vidéo](https://www.youtube.com/watch?v=-Xn1nmHjnHU) explique bien, étape par étape, ce que représente le test du $\chi^2$

```{r chi2}
tab <- table(dat$REG, dat$ZAU2)
chisq.test (tab)
```

La p-value est inférieure à 5%, on a donc une probabilité très faible de nous tromper en rejetant l'hypothèse d'indépendance. Donc on rejette l'indépendance, c'est à dire qu'on peut considérer que les variables sont liées.

Cependant, on a des effectifs nuls : il faut donc regrouper des modalités pour satisfaire aux conditions d'utilisation du test. Ici on écarte aussi les outre-mer (codes *01* à *04*) qui ont des effectifs faibles.

```{r}
dat_ss_om <- dat %>%
  filter (!REG %in% c('01', '02', '03', '04')) %>% 
  droplevels () %>% 
  mutate (ZAU3 = fct_recode (ZAU2, urbain = "111 - Grand pôle",
                             urbain = "211 - Moyen pôle", 
              urbain = "221 - Petit pôle", 
              "periurbain ou rural" = "112 - Couronne GP",
              "periurbain ou rural" = "212 - Couronne MP", 
              "periurbain ou rural" = "120 - Multipol grandes AU",
              "periurbain ou rural" = "300 - Autre multipol.", 
              "periurbain ou rural" = "222 - Couronne PP",
              "periurbain ou rural" = "400 - Commune isolée"))

levels (dat_ss_om$ZAU3)

tab <- table(dat_ss_om$REG, dat_ss_om$ZAU3)
chisq.test (tab)
```



### Un couple de variables qualitatives est-il "lié plus fortement" qu'un autre couple ?

- La statistique du $\chi^2$ dépend du nombre d'observations n et du nombre de modalités des variables. On ne peut donc pas comparer la force du lien entre plusieurs couples de variables
- Le **V de Cramer** quantifie l'intensité du lien entre deux variables qualitatives
- $V = \sqrt{\frac{\frac{\chi^2}{n}}{min(k-1,r-1)}} \in [0;1]$
- Analogie avec le coefficient de corrélation

```{r Cramer}
cramersV (tab)      # package lsr
```

On peut comparer cette grandeur pour plusieurs couples de variables qualitatives.

## Sur un croisement quantitatif/qualitatif

On cherche à comparer les moyennes d'une variable de deux groupes indépendants. La variable à comparer est quantitative, les deux (ou plus) groupes peuvent être vus comme issus d'une variable quantitative. Par exemple, on peut s'intéresser à la comparaison de **la densité (variable quantitative)** selon le **département (variable qualitative)**.

### 2 échantillons indépendants
 On ne regarde que 2 départements pour le moment. 

Exemple : On veut comparer la densité (log-transformée) des communes du Haut Rhin (68) et de celles du Val de Marne (94). Ce sont deux groupes indépendants (les communes du Haut Rhin ne sont pas liées à celles du Val de Marne).


#### Distribution normale
Si les deux groupes à observer sont normalement distribués, on se servira du *test de Student* (`t.test()` sur R).


```{r}
# travail préalable : on peut visualiser les données
data_test <- dat %>% filter(DEP %in% c(94, 68)) %>% 
  select(CODGEO, log_dens, densite, DEP)

ggplot(data = data_test, aes(x = DEP, y = log_dens, fill = DEP)) + 
  geom_boxplot() + theme (legend.position = "none")
```

```{r}
res <- t.test(log_dens ~ DEP, data = data_test, var.equal = TRUE)
res
```
La p-value est inférieure à 0.05. La moyenne de la log-densité dans le Val de Marne est significativement différente que celle dans le Haut Rhin.



#### Distribution non normale
Dans le cas où les distributions des groupes ne sont pas normales, on se servira du *test de Wilcoxon* (`wilcox.test()` sur R). Le test de Wilcoxon se base sur les rangs des valeurs des observations des deux groupes, ce qui permet de s'affranchir de l'hypothése de normalité.
Si on regarde directement la densité :

```{r}
res <- wilcox.test(densite ~ DEP, data = data_test)
res
```
La p-value est inférieure à 0.05. Les moyennes sont donc significativement différentes entre les deux départements.


### 2 échantillons dépendants
Cette fois, on cherche à comparer les moyennes entre deux groupes dépendants. On a alors 2 variables par observations. Par exemple, on veut comparer la population moyenne en 2014 par rapport à celle en 2009. Pour chaque commune, on dispose de la population en 2009 et 2014.
On se servira des mêmes tests que ceux vus précedemment, et on précisera l'argument `paired = TRUE`.


#### Distribution normale
Si les distributions sont normales, on se servira du *test de Student* 


```{r}
# travail préalable : on peut visualiser les données
data_test_d <- dat %>% filter(DEP == 94) %>% 
  select(CODGEO, P14_POP, P09_POP) %>% 
  gather(key = annee, value = population, -CODGEO) %>% 
  mutate(log_pop = log(population),
         annee = ifelse(annee == "P14_POP", "2014", "2009")) %>% 
  arrange(desc(annee), CODGEO)

g <- ggplot(data = data_test_d, aes(x = annee, y = log_pop, color = CODGEO)) + 
  geom_line(aes(group = CODGEO)) + geom_point() +
  theme (legend.position = "none")
ggplotly(g) # ajout de dynamisme 
```

```{r}
res <- t.test(log_pop ~ annee, data = data_test_d, paired = TRUE)
res
```

#### Distribution non normale
Si les distributions ne sont pas normales, on se servira du *test de Wilcoxon* 


```{r}
# travail préalable : on peut visualiser les données
g <- ggplot(data = data_test_d, aes(x = annee, y = population, color = CODGEO)) + 
  geom_line(aes(group = CODGEO)) + geom_point() +
  theme (legend.position = "none")
ggplotly(g)
```

```{r}
res <- wilcox.test(population ~ annee, data = data_test_d, paired = TRUE)
res

```
La p-value est inférieure à 0.05. La moyenne des populations est significativement différente d'une année sur l'autre. 



### Plusieurs échantillons indépendants
Si on cherche à comparer plus de 2 groupes, on doit se servir d'autres tests que ceux vus précedemment. De la même manière, on dispose d'un test qui fonctionne avec hypothése de normalité, et un qui permettra de s'affranchir de cette contrainte.


#### Echantillons normalement distribués (ANOVA)
Dans le cas où les données suivent une distribution normale, on se servira alors de *l'ANOVA* (Analysis of Variance). La fonction permettant d'effectuer une ANOVA sur R est `aov()`. Si on veut comparer la densité dans 3 départements différents, on fera : 

```{r anova_ggplot}
# travail préalable : on peut visualiser les données
data_test_m <- dat %>% filter(DEP %in% c(94, 68, 29)) %>% 
  select(CODGEO, densite, log_dens, DEP)

ggplot(data = data_test_m, aes(x = DEP, y = log_dens, fill = DEP)) +
  geom_boxplot() + theme (legend.position = "none")
```



```{r}
res <- aov(log_dens ~ DEP, data = data_test_m)
summary(res)
```

La p-value est inférieure à 0.05 : il y a bien une différence significative entre les groupes. Pour complèter l'analyse, il faudrait ensuite repartir sur les tests entre deux échantillons. Pour connaître les groupes deux à deux distincts, on peut aussi se servir de la fonction `TukeyHSD()`, qui calcule les comparaisons entre les groupes.


```{r}
TukeyHSD(res)
```
Les p-value étant toutes inférieures à 0.05, la différence est donc significative entre chaque groupes pris deux à deux.



Une autre représentation graphique utilisée pour une ANOVA est le graphique `errorbar`. 

```{r}
n <- nrow(data_test_m)

# creation des données pour chaque DEP
graph <- data_test_m %>% group_by(DEP) %>% 
  summarise(sd = sd(log_dens, na.rm = T),
            means = mean(log_dens, na.rm = T)) %>% 
  mutate(se = 1.96 * sd / sqrt(n))

ggplot(graph, aes(x = DEP, y = means, fill = DEP)) + 
  geom_bar(stat="identity", color="black", 
           position = position_dodge()) +
  geom_errorbar(aes(ymin = means - se, ymax = means + se), 
                width = 0.2) +
  theme (legend.position = "none")

```

#### Echantillons non normalement distribués (Kruskal-Wallis)
Si les distributions ne sont pas normales, on se servira du *test de Kruskal-Wallis*. Il généralise le test de Wilcoxon sur plus de deux groupes. Il fonctionne de manière similaire, en comparant les rangs des valeurs des différents groupes.

```{r kruskal_ggplot}
# travail préalable : on peut visualiser les données
ggplot(data = data_test_m, aes(x = DEP, y = densite, fill = DEP)) + geom_boxplot()
```

```{r}
res <- kruskal.test(densite ~ DEP, data = data_test_m)
res
```

La p-value est inférieure à 0.05 : il y a bien une différence significative entre les groupes. Pour complèter l'analyse, il faudrait ensuite repartir sur les tests entre deux échantillons. Pour connaître les groupes deux à deux distincts, on peut aussi se servir de la fonction `pariwise.wilcox.test`, qui calcule les comparaisons entre les groupes.

```{r}
pairwise.wilcox.test(data_test_m$densite, data_test_m$DEP,
                     p.adjust.method = "BH")

```

Ici, tous les départements sont significativement différents deux à deux, car les *p-values* sont toutes inférieures à 0,05.


## Résumé

Afin de savoir quel test s'applique dans chaque situation :

![](images/guide_tests.png)

```{r, echo= FALSE}

# fichiers dans Graphiques_tests
library(visNetwork)

point <- read.csv2(file = "Graphiques_tests/noeuds.csv", encoding = "UTF-8")
lien <- read.csv2(file = "Graphiques_tests/liens.csv", encoding = "UTF-8")


graph <- visNetwork(point, lien, height = '1000px', width = '100%') %>% 
  visNodes(shape = "box") %>%                        # square for all nodes
  visEdges(arrows = "to", smooth = F) %>%                            # arrow "to" for all edges
  visGroups(groupname = "question", color = "orange", shape = "ellipse") %>%   
  visGroups(groupname = "reponse", color = "lightgreen") %>%  
  visGroups(groupname = "test", color = "red") %>%    
  visGroups(groupname = "depart", color = "darkorchid", shape = "ellipse") %>%   
  visGroups(groupname = "interrogation", color = "lightslateblue") %>%    
  visOptions(highlightNearest = T,
             manipulation = F) %>%
  visHierarchicalLayout(direction = "LR", levelSeparation = 300) 


graph

```