10-Variance-I.Rmd

# Variance I {#variance}

```{r setup, include=FALSE, echo=FALSE, message=FALSE, results='hide'}
SciViews::R("infer", "model", lang = "fr")
```

<!-- Here, we need also to look at Type II versus type III + packages effectsize and afex. Discuss about eta^2, omega^2, Cohen's f etc (effect size in percent of variance) + the absolute effect. Also a power analysis with the pwr package and a way to estimate best n given a goal. Finally, one or two Shiny apps that allow to play with partition of variance to feel it better. A supplemental section that presents alternate homoscedasticity tests (Levene, ...). For effect size: https://cran.r-project.org/web/packages/effectsize/vignettes/anovaES.html -->

##### Objectifs {.unnumbered}

-   Pouvoir comparer les moyennes de plus de deux populations simultanément en utilisant des techniques de décomposition de la variance

-   Découvrir le modèle linéaire, anciennement analyse de variance (ANOVA)

-   Être capable d'effectuer des tests de comparaison multiples

-   Maîtriser l'équivalent non paramétrique à un facteur (test de Kruskal-Wallis)

##### Prérequis {.unnumbered}

Ce module continue la comparaison de moyennes entamée, pour deux populations, au module \@ref(moyenne). Assurez-vous d'avoir bien compris le test *t* de Student et les subtilités des tests d'hypothèse avant d'entamer la présente section.

## Le danger des tests multiples

Avant toute chose, nous préparons notre session R. Nous avons besoin de `SciViews::R` bien entendu, mais aussi des extensions `"infer"` et `"model"`. Tout cela est chargé en mémoire comme suit (n'oubliez pas ces codes dans un "chunk" de "setup" en début de vos documents Quarto ou R Markdown) :

```{r}
SciViews::R("infer", "model", lang = "fr")
```

Les tests *t* de Student et de Wilcoxon/Mann-Whitney sont limités à la comparaison de deux populations. Poursuivons notre analyse des crabes *L. variegatus*. Rappelez-vous, nous avons deux variétés (variable `species`, `B` pour bleue et `O`pour orange). Si nous voulons comparer simultanément les mâles et les femelles des deux variétés, cela nous fait quatre sous-populations à comparer (nous utilisons ici la fonction `paste()` qui rassemble des chaînes de caractère avec le tiret comme caractère séparateur `sep = "-"` pour former une variable facteur à quatre niveaux, `B-F`, `B-M`, `O-F`, `O-M`). Nous transformons cette variable en facteur à l'aide de `factor()` et nous ajoutons un label avec `labelise()`.

```{r}
crabs <- read("crabs", package = "MASS", lang = "fr")
crabs2 <- smutate(crabs, group  = labelise(
    factor(paste(species, sex, sep = "-")),
    "Groupe espèce - sexe", units = NA))
```

La Fig. \@ref(fig:crabs-rear) montre la largeur à l'arrière de la carapace chez les quatre groupes ainsi individualisés. Une représentation graphique adéquate avant de réaliser notre analyse lorsque le nombre de réplicats est important est le graphique en violon sur lequel nous superposons au moins les moyennes, et de préférence, les points également. Si le nombre de réplicats est plus faible, mais toujours supérieur à 7-8, nous pourrions utiliser le même type de graphique, mais avec des boites de dispersion plutôt (voir plus loin, Fig. \@ref(fig:anova2)). Avec encore moins de réplicats nous présenterons les points et les moyennes uniquement.

```{r crabs-rear, fig.cap="Largeur arrière en fonction du groupe de crabes *L. variegatus*. Graphique adéquat pour comparer les moyennes et distributions dans le cas d'un nombre important de réplicats (moyennes en rouge + observations individuelles en noir semi-transparent superposées à des graphiques en violon)."}
chart(data = crabs2, rear ~ group) +
  geom_violin() +
  geom_jitter(width = 0.05, alpha = 0.5) +
  stat_summary(geom = "point", fun = "mean", color = "red", size = 3)
```

##### À vous de jouer ! {.unnumbered}

`r h5p(227, height = 270, toc = "Graphique associé à l'ANOVA")`

Nous en profitons également pour essayer l'astuce proposée au module précédent. Au lieu de travailler sur la variable `rear` seule, nous allons étudier l'"aspect ratio" entre largeur à l'arrière (`rear`) et largeur maximale (`width`) de la carapace afin de nous débarrasser d'une source de variabilité triviale qui est qu'un grand crabe est grand pour toutes ses mensurations par rapport à un petit crabe. Nous prenons soin également de libeller cette nouvelle variable correctement avec `labelise()`. Enfin, nous ne conservons que les variables `species`, `sex`, `group` et `aspect` avec `sselect()` :

```{r}
crabs2 %>.%
  smutate(., aspect = labelise(
    as.numeric(rear / width),
    "Ratio largeur arrière / max", units = NA)) %>.%
  sselect(., species, sex, group, aspect) ->
  crabs2
skimr::skim(crabs2)
```

Nous avons 50 individus dans chacun des quatre groupes. **Lorsqu'il y a le même nombre de réplicats dans tous les groupes, on appelle cela un plan balancé**. C'est une situation optimale. Nous devons toujours essayer de nous en rapprocher le plus possible, car, si le nombre d'individus mesurés diffère fortement d'un groupe à l'autre, nous aurons forcément moins d'information disponible dans le ou les groupes moins nombreux, ce qui déforcera notre analyse.

##### À vous de jouer ! {.unnumbered}

`r h5p(228, height = 270, toc = "Plan balancé d'un expérience")`

Nous voyons également que la variable `aspect` semble avoir une distribution bimodale d'après le petit histogramme représenté dans le résumé. La Fig. \@ref(fig:anova0) avec `aspect` montre une différence plus marquée qu'avec `rear` seul (voir Fig. \@ref(fig:crabs-rear)).

```{r anova0, fig.cap="Ratio largeur arrière/largeur max en fonction du groupe de crabes *L. variegatus*. Graphique adéquat pour comparer les moyennes et distributions dans le cas d'un nombre important de réplicats (moyennes en rouge + observations individuelles en noir semi-transparent superposées à des graphiques en violon)."}
chart(data = crabs2, aspect ~ group) +
  geom_violin() +
  geom_jitter(width = 0.05, alpha = 0.5) +
  stat_summary(geom = "point", fun = "mean", color = "red", size = 3)
```

Nous voyons ici beaucoup mieux que la distribution bimodale qui avait été observée sur l'ensemble des données est essentiellement due au dimorphisme sexuel plutôt qu'à des différences entre les variétés. Mais qu'en est-il plus précisément ? Si nous regardons attentivement, il semble que les moyennes pour les crabes bleus soient légèrement inférieures à celles des crabes oranges.

Comment comparer valablement ces quatre groupes ? Comme nous savons maintenant comparer deux groupes à l'aide d'un test *t* de Student, il est tentant d'effectuer toutes les comparaisons deux à deux et de résumer l'ensemble, par exemple, dans un tableau synthétique. Cela fait quand même beaucoup de comparaisons (`B-F` \<-\> `B-M`, `B-F` \<-\> `O-F`, `B-F` \<-\> `O-M`, `B-M` \<-\> `O-F`, `B-M`\<-\> `O-M`, et finalement `O-F` \<-\> `O-M`). Cela fait donc six comparaisons à réaliser.

N'oublions pas que, à chaque test, nous prenons un risque de nous tromper. **Le risque de se tromper au moins une fois dans l'ensemble des tests est alors décuplé en cas de tests multiples.** Prenons un point de vue naïf, mais qui suffira ici pour démontrer le problème qui apparaît. Admettons que le risque de nous tromper est constant, que nous rejetons ou non $H_0$ (donc les risques $\alpha$ et $\beta$ sont identiques), et qu'il est de l'ordre de 10% dans chaque test individuellement[^10-variance-i-1]. La seule solution acceptable est que les interprétations de *tous* les tests soient correctes. Considérant chaque interprétation indépendante, nous pouvons multiplier les probabilités d'avoir un test correct (90%) le nombre de fois que nous faisons le test, soit $0,9 \times 0,9 \times 0,9 \times 0,9 \times 0,9 \times 0,9 = 0,9^6 = 0,53$. Tous les autres cas ayant au moins un test faux, nous constatons que notre analyse globale sera incorrecte $1 - 0,53 = 47\%$ du temps[^10-variance-i-2]. **Notre analyse sera incorrecte une fois sur deux environ.**

[^10-variance-i-1]: Attention ! vous savez bien que c'est plus compliqué que cela. D'une part, le risque de se tromper est probablement différent si on rejette $H_0$ ($\alpha$) ou non ($\beta$), et ces risques sont encore à moduler en fonction de la probabilité *a priori*, un cas similaire au dépistage d'une maladie plus ou moins rare, rappelez-vous, au module \@ref(proba).

[^10-variance-i-2]: Dans R, vous pouvez utiliser `choose(n, j)` pour calculer le coefficient binomial. Donc votre calcul du risque de se tromper au moins une fois dans un ensemble de `n` tests dont le risque individuel est `r` sera `1 - (1 - r)^choose(n, 2)`.

```{block, type='info'}
De manière générale, le nombre de combinaisons deux à deux possibles dans un set de `n` groupes distincts sera calculé à l'aide du coefficient binomial que nous avions déjà rencontré avec la distribution du même nom, ici avec $j$ valant deux.

$$C^j_n = \frac{n!}{j!(n-j)!}$$

Toujours avec notre approche naïve du risque d'erreur individuel pour un test $r$ de 10%, le risque de se tromper au moins une fois est alors :

$$1 - (1 - r)^{C^2_n}$$

Voici ce que cela donne comme risque de se tromper dans au moins un des tests en fonction du nombre de groupes à comparer :

|Groupes comparés 2 à 2  |  2  |  3  |  4  |  6  |  8  |  10 |
|:-----------------------|:---:|:---:|:---:|:---:|:---:|:---:|
|Risque individuel = 10% | 10% | 27% | 47% | 79% | 95% | 99% |

Clairement, on oublie cette façon de faire ! Prendre le risque de se tromper 99 fois sur 100 en comparant 10 groupes différents n'est pas du tout intéressante comme perspective.

```

Nous allons donc travailler différemment... Ci-après nous verrons qu'une simplification des hypothèses et l'approche par décomposition de la variance est une option bien plus intéressante (ANalysis Of VAriance ou ANOVA). Ensuite, nous reviendrons vers ces comparaisons multiples deux à deux, mais en prenant des précautions pour éviter l'inflation du risque global de nous tromper.

##### À vous de jouer ! {.unnumbered}

`r h5p(229, height = 270, toc = "Comparaison deux à deux")`

## ANOVA à un facteur

Au lieu de nous attaquer aux comparaisons deux à deux, nous pouvons considérer une hypothèse unique que les moyennes de $k$ populations (nos quatre groupes différents de crabes, par exemple) sont toutes égales. L'hypothèse alternative sera qu'au moins une des moyennes diffère des autres (sans préciser laquelle). En formulation mathématique, cela donne :

-   $H_0: \mu_1 = \mu_2 = ... = \mu_k$

-   $H_1: \exists(i, j) \mathrm{\ tel\ que\ } \mu_i \neq \mu_j$

Notre hypothèse nulle est très restrictive, mais par contre, l'hypothèse alternative est très vague, car nous ne savons pas **où** sont les différences à ce stade si nous rejetons $H_0$. Nous nous en occuperons plus tard.

```{block, type='info'}
**Propriété d'additivité des parts de variance**. La variance se calcule comme :

$$var_x = \frac{SCT}{ddl}$$

Avec $SCT$, la somme des carrés totaux, soit $\sum_{i = 1}^n (x_i - \bar{x})^2$, la somme des carrés des écarts à la moyenne générale. Les ddl sont les **degrés de liberté** déjà rencontrés à plusieurs reprises qui valent $n - 1$ dans le cas de la variance d'un échantillon.

Cette variance peut être *partitionnée*. C'est-à-dire que, si la variance totale se mesure d'un point A à un point C, l'on peut mesurer la part de variance d'un point A à un point B, puis d'autre part, d'un point B à un point C, et dans ce cas,

$$SCT = SC_{A-C} = SC_{A-B} + SC_{B-C}$$

Cette propriété, dite d'additivité des variances, permet de décomposer la variance totale à souhait tout en sachant que la somme des différentes composantes donne toujours la même valeur que la variance totale.

```

##### À vous de jouer ! {.unnumbered}

`r h5p(100, height = 400, toc = "Hypothèse nulle de l'ANOVA")`

### Modèle de l'ANOVA

Mais qu'est-ce que cette propriété d'additivité des variances vient faire ici ? Nous souhaitons comparer des moyennes, non ? Effectivement, mais considérons le **modèle mathématique suivant :**

$$y_{ij} = \mu + \tau_j + \epsilon_i \mathrm{\ avec\ } \epsilon \sim N(0, \sigma^2)$$

Avec l'indice $j = 1 .. k$ populations et l'indice $i = 1 .. n$ observations du jeu de données. Chaque observation $y_{ij}$ correspond à deux écarts successifs de la moyenne globale $\mu$ : une constante "tau" par population $\tau_j$ d'une part et un terme $\epsilon_i$ que l'on appelle les **résidus** et qui est propre à chaque observation individuelle. C'est ce dernier terme qui représente la partie statistique du modèle avec une distribution normale centrée sur zéro et avec une variance $\sigma^2$ que nous admettrons constante et identique pour toutes les populations par construction.

Le graphique à la Fig. \@ref(fig:anova1) représente une situation typique à trois sous-populations.

```{r anova1, echo = FALSE, fig.cap="Décomposition de la variance dans un cas à trois populations A, B et C fictives."}
dat <- dtx(Population = c("A", "A", "A", "B", "B", "B", "C", "C", "C"),
  Réponse = c(2.2, 2.4, 3.3, 4.9, 5.2, 6.1, 4.3, 6.4, 6.5))
means <- c(mean(dat$Réponse), tapply(dat$Réponse, dat$Population, mean))
annot <- dtx(x = c(1.1, 1.4, 1.7, 2.1, 2.4, 2.7),
  ystart = c(means[1], means[1], means[2], means[1], means[1], means[3]),
  yend = c(2.2, means[2], 2.2, 4.9, means[3], 4.9),
  color = c(2L, 2L, 2L, 3L, 3L, 3L))
annot2 <- dtx(xstart = c(0.95, 1.95), xend = c(1.75, 2.75),
  y = c(2.2, 4.9), color = 2:3)
library(ggplot2)
chart(data = dat, Réponse ~ Population %col=% Population) +
  geom_point() +
  theme(legend.position = "none") +
  geom_hline(yintercept = means, col = 1:4, linetype = "dashed", size = 0.5) +
  geom_segment(data = annot2, aes(x = xstart, xend = xend, y = y , yend = y),
    color = annot2$color, linetype = "dotted", size = 0.5) +
  geom_segment(data = annot, aes(x = x, xend = x, y = ystart, yend = yend),
    color = annot$color, size = 0.5, arrow = arrow(length = unit(0.2, "cm"))) +
  #annotate("text", x = 1, y = means[1] + 0.2, color = 1,
  #  label = "moyenne ~ générale ~ mu", parse = TRUE) +
  annotate("text", x = 1, y = means[1] + 0.2, color = 1,
    label = "moyenne ~ globale ~ mu", parse = TRUE) +
  annotate("text", x = 3.1, y = means[2] + 0.2, color = 2,
    label = "moyenne pour A") +
  annotate("text", x = 3.1, y = means[3] + 0.2, color = 3,
    label = "moyenne pour B") +
  annotate("text", x = 3.1, y = means[4] + 0.2, color = 4,
    label = "moyenne pour C") +
  annotate("text", x = 1.25, y = 3.25, color = 2, label = "total") +
  annotate("text", x = 1.55, y = 3.6, color = 2, label = "inter ~ tau[A]", parse = TRUE) +
  annotate("text", x = 1.85, y = 2.5, color = 2, label = "intra ~ epsilon[i]", parse = TRUE) +
  annotate("text", x = 2.25, y = 4.8, color = 3, label = "total") +
  annotate("text", x = 2.55, y = 5.1, color = 3, label = "inter ~ tau[B]", parse = TRUE) +
  annotate("text", x = 2.85, y = 5.2, color = 3, label = "intra ~epsilon[i]", parse = TRUE)
```

Notons que ce modèle à trois termes représente bien la situation qui nous intéresse, mais aussi, qu'il décompose la variance totale (entre $\mu$ et chaque point observé) en deux : ce que nous appellerons le terme **inter** représentant l'écart entre la moyenne globale $\mu$ et la moyenne de la sous-population concernée ($\tau_j$) et le terme **intra** depuis cette moyenne de la sous-population jusqu'au point observé ($\epsilon_i$).

##### À vous de jouer ! {.unnumbered}

`r h5p(99, height = 400, toc = "Variance intergroupe dans l'ANOVA")`

D'une part, nous nous trouvons dans une situation d'additivité de la variance si nous décidons de calculer ces "variance inter" et "variance intra". D'autre part, sous $H_0$ nous sommes sensés avoir toutes les moyennes égales à $\mu$, et donc, tous les $\tau_j = 0$. Donc, les valeurs non nulles de $\tau_j$ ne doivent qu'être dues au hasard de l'échantillonnage et être par conséquent largement inférieures à la variabilité entre les individus, ou variance intra $\epsilon_i$. La Fig. \@ref(fig:anova2) représente deux cas avec à gauche une situation où $H_0$ est plausible, et à droite une situation où elle est très peu plausible. Notez qu'à gauche la variation entre les observations (intra) est bien plus grande que l'écart entre les moyennes (inter), alors qu'à droite c'est l'inverse.

```{r anova2, echo=FALSE, fig.cap="A. Cas fictif avec moyennes probablement égales entre populations (étalement des points bien plus large que l'écart entre les moyennes), B. cas où les moyennes sont probablement différentes (écart des moyennes \"inter\" bien plus grand que l'étalement des points en \"intra\"-population)."}
set.seed(25431)
dat1 <- dtx(
  Y = rnorm(30, sd = 0.3),
  Population = rep(c("A", "B", "C"), each = 10))
dat2 <- dtx(
  Y = rnorm(30, sd = 0.1) + rep(c(0, -1, 1), each = 10),
  Population = rep(c("A", "B", "C"), each = 10))
pl1 <- chart(data = dat1, Y ~ Population) +
  geom_boxplot() +
  geom_jitter(width = 0.05) +
  stat_summary(geom = "point", fun = "mean", color = "red", size = 3) +
  ylim(-1.5, 1.5)
pl2 <- chart(data = dat2, Y ~ Population) +
  geom_boxplot() +
  geom_jitter(width = 0.05) +
  stat_summary(geom = "point", fun = "mean", color = "red", size = 3) +
  ylim(-1.5, 1.5)
combine_charts(list(pl1, pl2))
```

Intuitivement, une comparaison de "inter" et "intra" permet de différencier la situation de gauche de celle de droite dans la Fig. \@ref(fig:anova2). Si cette comparaison est faite sous forme d'un ratio "inter"/"intra", alors ce ratio sera faible et tendra vers zéro sous $H_0$ (cas A), alors qu'il sera d'autant plus élevé que $H_0$ devient de moins en moins plausible (cas B).

##### Calcul de l'ANOVA {.unnumbered}

Calcul des sommes des carrés (inter- et intragroupes). Considérant :

-   *i* = indice des observations au sein du jeu de données de 1 à *n*,
-   *j* = facteurs (sous-populations de 1 à *k*),
-   $\bar{y}$ = moyenne générale de l'échantillon,
-   $\bar{y_j}$ = moyenne de la *j*^ème^ population.

La somme des carrés inter $SC_{inter}$ et la somme des carrés intra $SC_{intra}$ se calculent comme suit :

$$
\begin{aligned}
SC_{inter} = \sum_{i=1}^n{(\bar{y_j} - \bar{y})^2} && SC_{intra} = \sum_{i=1}^n{(y_{ij} - \bar{y_j})^2}
\end{aligned}
$$

À ces sommes des carrés, nous pouvons associer les degrés de liberté suivants :

-   *k* -- 1 pour l'intergroupe
-   *n* -- *k* pour l'intragroupe

Sachant que les parts de variances sont les $\frac{SC}{ddl}$ et sont appelés "carrés moyens", nous construisons ce qu'on appelle le **tableau de l'ANOVA** de la façon suivante :

| Type            |    Ddl    | Somme carrés |   Carré moyen (*CM*)   | Statistique *F*~obs~  | P (\>*F*) |
|:----------------|:---------:|:------------:|:----------------------:|:---------------------:|:---------:|
| Inter (facteur) |  *k* - 1  | *SC~inter~*  | *SC~inter~/ddl~inter~* | *CM~inter~/CM~intra~* |    ...    |
| Intra (résidus) | *n* - *k* | *SC~intra~*  | *SC~intra~/ddl~intra~* |                       |           |

La **statistique *F*~obs~ est le rapport des carrés moyens inter/intra.** Elle représente donc le ratio que nous avons évoqué plus haut comme moyen de quantifier l'écart par rapport à $H_0$. Le test consiste à calculer la valeur *P* associée à cette statistique. Pour cela, il nous faut une distribution statistique théorique de *F* sous $H_0$. C'est un biologiste - statisticien célèbre nommé Ronald Aylmer Fisher qui l'a calculée. Pour cette raison, elle s'appelle la distribution *F* de Fisher.

##### À vous de jouer ! {.unnumbered}

`r h5p(101, height = 400, toc = "Analyse du tableau de l'ANOVA")`

### Distribution *F*

La distribution *F* est une distribution asymétrique n'admettant que des valeurs nulles ou positives, d'une allure assez similaire à la distribution du $\chi^2$ que nous avons étudié au module \@ref(chi2). Elle est appelée loi de Fisher, ou encore, loi de Fisher-Snedecor. Elle a une asymptote horizontale à $+\infty$. La distribution *F* admet deux paramètres, respectivement les degrés de liberté au numérateur (inter) et au dénominateur (intra). La Fig. \@ref(fig:fplot) représente la densité de probabilité d'une loi *F* typique[^10-variance-i-3].

[^10-variance-i-3]: Les fonctions qui permettent les calculs relatifs à la distribution *F* dans R sont `<q|p|r|d>f()`. En outre, `dist_f()` génère un objet permettant d'en tracer facilement le graphique. Leur utilisation est similaire à celle des distributions vues au module \@ref(proba).

<!-- et les snippets correspondants dans la SciViews Box sont disponibles à partir de `.if`. -->

```{r fplot, echo=FALSE, warning=FALSE, fig.cap = "Allure typique de la densité de probabilité de la distribution F (ici ddl inter = 5 et ddl intra = 20). Plus *F~obs* est grand, plus l'hypothèse nulle est suspecte. La zone de rejet est donc positionnée à droite (en rouge)."}
# Fisher-Snedecor's F distribution (density probability) with parameter:
#.df1 <- 5; .df2 <- 20 # numerator (.df1) and denominator (.df2) df
#.col <- 1; .add <- FALSE # Plot parameters
#.x <- seq(0, qf(0.999, df1 = .df1, df2 = .df2), l = 1000)  # Quantiles
#.d <- function(x) df(x, df1 = .df1, df2 = .df2)            # Distribution function
#.q <- function(p) qf(p, df1 = .df1, df2 = .df2)            # Quantile for lower-tail prob
#.label <- bquote(F(.(.df1), .(.df2)))                      # Curve parameters
##curve(.d(x), xlim = range(.x), xaxs = "i", n = 1000, col = .col,
##  add = .add, xlab = "Quantiles", ylab = "Probability density") # Curve
##abline(h = 0, col = "gray") # Baseline

#q_ref <- qf(0.05, df1 = .df1, df2 = .df2, lower.tail = FALSE)
#.x2 <- .x
#.x2[.x2 < q_ref] <- NA

#chart(data = dtx(Quantiles = .x, Prob = .d(.x)), Prob ~ Quantiles) +
#  geom_hline(yintercept = 0, col = "gray") +
#  geom_ribbon(aes(x = .x, ymin = 0, ymax = .d(.x)), fill = "gray", alpha = 0.2) +
#  geom_ribbon(aes(x = .x2, ymin = 0, ymax = .d(.x2)), fill = "red", alpha = 0.2) +
#  geom_line() +
#  xlab("Quantile") +
#  ylab("Densité de probabilité") +
#  geom_vline(xintercept = 3.5, col = "red") +
#  annotate("text", x = 1.2, y = 0.1, label = "Zone de non rejet", col = "black") +
#  annotate("text", x = 5, y = 0.1, label = "plain(Zone ~ de ~ rejet) == plain(seuil) ~ alpha", parse = TRUE, col = "red") +
#  annotate("segment", x = 4.5, y = 0.09, xend = 4.2, yend = 0.01,
#    size = 0.5, col = "red", arrow = arrow(length = unit(0.2, "cm"))) +
#  annotate("text", x = 3.8, y = 0.35, label = "F[obs]", parse = TRUE,
#    col = "red", fontface = "bold.italic") +
#  annotate("text", x = 4.7, y = 0.3, label = "(aire à droite = valeur P)",
#    col = "red", fontface = "bold.italic")

f <- dist_f(df1 = 5, df2 = 20)
chart(f) +
  geom_funfill(fun = dfun(f), from = quantile(f, 1 - 0.05), to = 8) +
  geom_hline(yintercept = 0, col = "gray") +
  geom_vline(xintercept = 3.5, col = "red") +
  ylab("Densité de probabilité") +
  xlim(-0.1, 7) +
  annotate("text", x = 1.2, y = 0.1, label = "Zone de non rejet", col = "black") +
  annotate("text", x = 5, y = 0.1, label = "plain(Zone ~ de ~ rejet) == plain(seuil) ~ alpha", parse = TRUE, col = "red") +
  annotate("segment", x = 4.5, y = 0.09, xend = 4.2, yend = 0.01,
    size = 0.5, col = "red", arrow = arrow(length = unit(0.2, "cm"))) +
  annotate("text", x = 3.8, y = 0.35, label = "F[obs]", parse = TRUE,
    col = "red", fontface = "bold.italic") +
  annotate("text", x = 4.7, y = 0.3, label = "(aire à droite = valeur P)",
    col = "red", fontface = "bold.italic")
```

Nous commençons à avoir l'habitude maintenant. La valeur *P* est calculée comme l'aire à droite du quantile correspondant à *F~obs~*. Enfin, nous rejetterons $H_0$ seulement si la valeur *P* est inférieure au seuil $\alpha$ qui a été choisi préalablement au test. Ceci revient à constater que, graphiquement, *F~obs~* vient se positionner dans la zone de rejet en rouge comme sur la Fig. \@ref(fig:fplot).

```{block, type='info'}
Lorsque nous réalisons une ANOVA, nous travaillons avec deux variables. La variable qualitative sert à découper l'échantillon en plusieurs sous-populations, alors que la variable quantitative est utilisée pour calculer les moyennes. Leurs rôles diffèrent donc, et elles ne sont bien évidemment pas interchangeables. En statistique, nous parlerons de **variable réponse** pour la variable qui mesure l'effet étudié et **variable explicative** pour la variable qui permet d'expliquer la réponse observée (et dans le cas de l'ANOVA, c'est le découpage en sous-populations).

Si le modèle est représenté par une formule dans R, la variable réponse sera à la gauche de la formule, et la (ou les) variable(s) explicative(s) à droite : `var_réponse ~ var_explicative`.

```

##### Conditions d'application {.unnumbered}

-   échantillon représentatif (par exemple, aléatoire),
-   observations indépendantes,
-   variable dite **réponse** quantitative,
-   une variable dite **explicative** qualitative à trois niveaux ou plus,
-   distribution **normale** des résidus $\epsilon_i$,
-   **homoscédasticité** (même variance intragroupes, *"homoscedasticity"* en anglais, opposé à hétéroscédasticité = variance différente entre les groupes).

Les deux dernières conditions d'applications doivent être vérifiées. La normalité des résidus doit être rencontrée aussi bien que possible. Un graphique quantile-quantile des résidus permet de se faire une idée, comme sur la Fig. \@ref(fig:resid1). Néanmoins, le test étant relativement robuste à des petites variations par rapport à la distribution normale, surtout si ces variations sont symétriques, nous ne serons pas excessivement stricts ici.

```{r resid1, echo=FALSE, results='hide', fig.cap="Graphique quantile-quantile appliqué aux résidus d'une ANOVA pour déterminer si leur distribution se rapproche d'un loi normale."}
# Just a random dataset...
trees <- read("trees", package = "datasets")
ANOVA <- lm(data = trees, volume ~ diameter)
#anova(.anova)
#plot(.anova, which = 2)
#qq_slope <- function(sample) {
#  probs <- c(0.25, 0.75)
#  x <- qnorm(probs)
#  y <- quantile(sample, probs, na.rm = TRUE)
#  diff(y) / diff(x)
#}
#qq_intercept <- function(sample) {
#  probs <- c(0.25, 0.75)
#  x <- qnorm(probs)
#  y <- quantile(sample, probs, na.rm = TRUE)
#  slope <- diff(y) / diff(x)
#  y[1L] - slope * x[1L]
#}
#.anova %>.%
#  broom::augment(.) %>.%
#  chart(data = ., aes(sample = .std.resid)) +
#  geom_qq() +
#  #geom_qq_line() +
#  geom_abline(slope = qq_slope(.$.std.resid), intercept = qq_intercept(.$.std.resid),
#    colour = "darkgray") +
#  labs(x = "Theoretical quantiles", y = "Standardized residuals") +
#  ggtitle("Normal Q-Q")

#trees_anova2 <- broom::augment(trees_anova)
#invisible(car::qqPlot(.anova2[[".std.resid"]], distribution = "norm",
#  envelope = 0.95, col = "Black", xlab = "Quantiles théoriques (distri. normale)",
#  ylab = "Résidus"))
chart$qqplot(ANOVA, lang = "fr")
```

Pour rappel, ce graphique s'interprète comme suit : si les points s'alignent plutôt bien sur la droite, nous pouvons considérer que la distribution des résidus (ici la version standardisée, donc avec moyenne nulle et variance de un) est compatible avec la distribution théorique comparée qui est ici la distribution normale.

La condition d'homoscédasticité est plus sensible. Elle mérite donc d'être vérifiée systématiquement et précisément. Différents tests d'hypothèse existent pour le vérifier, comme le test de Bartlett, le test de Levene, etc. Nous vous proposons ici d'utiliser le test de Bartlett. Ses hypothèses sont :

-   $H_0: var_1 = var_2 = ... = var_k$ (homoscédasticité)
-   $H_1: \exists(i, j) \mathrm{\ tel\ que\ } var_i \neq var_j$ (hétéroscédasticité)

Si la valeur *P* est inférieure au seuil $\alpha$ fixé au préalable, nous devrons rechercher une transformation des variables qui stabilisera la variance. La première transformation à essayer en biologie et la transformation logarithmique surtout si les valeurs négatives de la variable réponse ne sont pas possibles, signe d'une distribution qui est peut être plutôt de type log-normale pour cette variable. Si aucune transformation ne stabilise la variance, nous devrons nous rabattre vers un test non paramétrique équivalent, le test de Kruskal-Wallis que nous aborderons plus loin dans ce module. Il n'existe pas d'équivalent de l'ANOVA avec variances inégales, comme le test de Welch *vs* test *t* de Student classique dans le cas de la comparaison de deux moyennes.

<!-- La SciViews Box propose des snippets pour accéder à ces différentes analyses. Dans le menu `hypothesis tests: variances` ou `.hv` nous trouvons trois tests dont celui de Bartlett. Dans le menu `hypothesis tests: means` ou `.hm` se trouvent les templates pour l'ANOVA, ainsi que les graphiques d'analyse des résidus dont le graphique quantile-quantile.-->

##### À vous de jouer ! {.unnumbered}

`r h5p(102, height = 400, toc = "Variance intragroupe dans l'ANOVA")`

##### Résolution de notre exemple {.unnumbered}

Nous commençons par déterminer si nous avons une homoscédasticité. Considérons un seuil $\alpha$ de 5% pour tous nos tests. Ensuite, nous utilisons `bartlett.test()` avec une formule qui représente bien la variable réponse *en fonction de* (le tilde `~`) la variable explicative, soit ici, `aspect ~ group` :

```{r}
bartlett.test(data = crabs2, aspect ~ group)
```

Nous rejetons $H_0$. Il n'y a pas homoscédasticité. Calculons par exemple le logarithme népérien de `aspect` et réessayons :

```{r}
crabs2 <- smutate(crabs2, log_aspect = ln(aspect))
bartlett.test(data = crabs2, log_aspect ~ group)
```

Ici cela ne fonctionne pas. Cela fait pire qu'avant. La transformation inverse (`exp()`) peut être essayée, mais ne stabilise pas suffisamment la variance. Après divers essais, il s'avère qu'une transformation puissance cinq stabilise bien la variance.

```{r}
crabs2 <- smutate(crabs2, aspect5 = aspect^5)
bartlett.test(data = crabs2, aspect5 ~ group)
```

La Fig. \@ref(fig:crabs1) montre la distribution dans les différents groupes de la variable transformée.

```{r crabs1, fig.cap="Transformation puissance cinq du ratio largeur arrière/largeur max en fonction du groupe de crabes *L. variegatus*."}
chart(data = crabs2, aspect5 ~ group) +
  geom_violin() +
  geom_jitter(width = 0.05, alpha = 0.5) +
  stat_summary(geom = "point", fun = "mean", color = "red", size = 3) +
  ylab("(ratio largeur arrière/max)^5")
```

Nous poursuivons sur une description des données utile pour l'ANOVA :

<!-- Un snippet dédié est disponible dans le menu `hypothesis tests: means` à partir de `.hm`. -->

```{r}
crabs2 %>.%
  sgroup_by(., group) %>.%
  ssummarise(.,
    mean  = fmean(aspect5),
    sd    = fsd(aspect5),
    count = fsum(!is.na(aspect5)))
```

Ensuite l'ANOVA proprement dite se réalise en deux étapes : (1) `lm()` calcule le modèle sur base d'une formule de type `var_réponse ~ var_explicative` (ici, `aspect5 ~ group`) et nous récupérons le résultat de ce calcul dans `crabs2_anova` pour utilisation ultérieure, et (2) la fonction `anova()` appliquée sur cet objet nous produit le tableau de l'ANOVA. Pour une version plus propre de ce tableau, nous utilisons enfin `tabularise()` :

```{r}
anova(crabs2_anova <- lm(data = crabs2, aspect5 ~ group)) |> tabularise()
```

Nous observons dans le tableau de l'ANOVA que la valeur *P* est très faible et inférieure à $\alpha$. Nous rejetons $H_0$. Nous pouvons dire que le ratio largeur arrière / max à la puissance cinq diffère significativement entre les groupes (ANOVA, F = 150, ddl = 3 & 196, valeur *P* \<\< 10^-3^).

Ce n'est qu'une fois l'ANOVA réalisée que nous pouvons calculer et visionner la distribution des résidus à partir de notre objet `crabs2_anova` (Fig. \@ref(fig:anova-resid)).

```{r anova-resid, echo=TRUE, fig.cap="Graphique quantile-quantile des résidus pour l'analyse de `aspect^5`."}
chart$qqplot(crabs2_anova, lang = "fr")
```

Cette distribution s'éloigne un peu aux extrêmes de la distribution normale. Nous sommes dans un cas un peu limite. Soit nous considérons que l'écart est suffisamment léger pour ne pas trop influer sur l'ANOVA, soit nous devons passer à un test non paramétrique. Nous ferons les deux ici à des fins de comparaison.

##### Pour en savoir plus {.unnumbered}

-   [L'ANOVA expliquée en trois minutes](https://www.youtube.com/watch?v=lITNHx2z5FE)

-   Introduction to ANOVA (en anglais). [Part I](https://youtu.be/QUQ6YppWCeg), [part II](https://youtu.be/fFnOD7KBSbw), [part III](https://youtu.be/XdZ7BRqznSA), [part IV](https://youtu.be/WUoVftXvjiQ), and [part V](https://youtu.be/kO8t_q-AXHE).

-   Explication de l'analyse de variance en détaillant le calcul par la Kahn Academy. [Partie I](https://www.youtube.com/watch?v=tjolTrwJhjM), [partie II](https://youtu.be/DMo9yofC5C8) et [partie III](https://youtu.be/y8nRhsixBPs). Assez long : près de 3/4h en tout. Ne regardez que si vous n'avez pas compris ce que sont les sommes des carrés.

## Test "post-hoc"

Avec l'ANOVA, lorsque nous rejetons $H_0$ comme dans le cas de notre exemple, nous savons qu'il y a au moins deux moyennes qui diffèrent l'une de l'autre, mais nous ne savons toujours pas lesquelles à ce stade. Notre analyse n'est pas terminée. Nous allons revisiter les tests de comparaison multiples deux à deux, mais en prenant des précautions particulières pour éviter l'inflation du risque d'erreur.

Tout d'abord, nous mettons en place un garde-fou. Nous effectuons *toujours* une ANOVA en premier lieu, et nous n'envisageons les comparaisons multiples que lorsque $H_0$ est rejetée. Cela évite beaucoup de faux positifs. On appelle cela des tests "post hoc" car ils ne sont pas planifiés d'emblée, mais suivent un conclusion préalable (ici, le rejet de $H_0$ dans un test ANOVA).

##### À vous de jouer ! {.unnumbered}

`r h5p(230, height = 400, toc = "tests post hoc")`

Une approche simple consisterait à modifier notre seuil $\alpha$ pour chaque test individuel vers le bas afin que le risque de se tromper dans au moins un des tests ne dépasse pas la valeur de $/alpha$ que nous nous sommes fixée. C'est la **correction de Bonferroni**. Elle consiste à diviser la valeur de $/alpha$ par le nombre de tests simultanés nécessaires, et d'utiliser cette valeur corrigée comme seuil $\alpha$ de chaque test de Student individuel. Dans le cas de quatre populations, nous avons vu qu'il y a six comparaisons multiples deux à deux. Donc, en appliquant un seuil corrigé de $\alpha/6 = 0,05 / 6 = 0,00833$ pour chaque test, on aura la probabilité suivante pour $1 - \alpha$ :

$$(1 - 0,05/6)^6 = 0,951$$

Donc, le risque *global* pour l'ensemble des six tests est bien de 1 - 0,951 = 0,049, soit 5%. Si elle a le mérite d'être simple, cette façon de faire n'est pas considérée comme la plus efficace. Actuellement, la méthode **HSD de Tukey** est préférée. HSD veut dire "Honest Significant Difference". La technique consiste à calculer l'écart minimal des moyennes pour considérer qu'elles sont significativement différentes l'une de l'autre. Ensuite, pour chaque comparaison deux à deux, l'écart entre les moyennes est comparé à cette valeur de référence[^10-variance-i-4]. Dans R, le test peut se visualiser soit sous une forme textuelle, soit sous une version graphique qui résume les comparaisons. Le code est un peu plus complexe et fait intervenir le package {multcomp}. Inspirez-vous de ce code dans vos projets.

[^10-variance-i-4]: Nous ne détaillons pas le calcul du test de Tukey dans le cours, mais vous pouvez aller voir [ici](https://www.itl.nist.gov/div898/handbook/prc/section4/prc471.htm).

<!-- Le test de comparaisons multiples de Tukey est accessible à partir des "snippets" dans la SciViews Box (`.hm` pour le menu `hypothesis tests: means`, puis `anova - multiple comparisons`). -->

```{r}
# Version textuelle avec summary()
summary(crabs2_posthoc <- confint(multcomp::glht(crabs2_anova,
  linfct = multcomp::mcp(group = "Tukey"))))
# Version graphique avec plot()
oma <- par(oma = c(0, 5.1, 0, 0))
plot(crabs2_posthoc)
par(oma)
rm(oma)
```

Avec la version textuelle, nous comparons comme à notre habitude les valeurs *P* renseignées dans la dernière colonne du tableau (nommée `Pr(>|t|)`) au seuil $\alpha$ avec la règle habituelle de rejet si *P* \< $\alpha$ (donc, différence des moyennes de cette paire significative au seuil $\alpha$ choisi). La comparaison se fait pour chaque ligne du tableau.

Avec la version graphique, les barres d'erreurs horizontales représentent les intervalles de confiance à $1 - \alpha$ autour de la différence des moyennes pour la paire considérée (indiquée en libellé de l'axe des ordonnées). Si cet intervalle de confiance ne contient pas zéro (par ailleurs mis en évidence par un trait vertical pointillé sur le graphique), nous pourrons considérer que les deux moyennes diffèrent de manière significative au seuil $\alpha$ choisi. Ici aussi, l'interprétation se fait ligne après ligne, c'est à dire, paire après paire.

Si nous gardons notre seuil $\alpha$ de 5%, seuls les mâles ne diffèrent (tout juste) pas entre eux (ligne `O-M - B-M` dans la version textuelle et dans le graphique) avec une valeur *P* de 6%. Toutes les autres comparaisons deux à deux sont significativement différentes.

En conclusion, tous les groupes diffèrent de manière significative sauf les mâles (test HSD de Tukey, valeur *P* \< 5%).

##### Conditions d'application {.unnumbered}

Les conditions d'application pour le test post hoc de Tukey sont les mêmes que pour l'ANOVA.

##### À vous de jouer ! {.unnumbered}

`r learnr("A10La_anova", title = "ANOVA et tests post-hoc", toc = "Analyse de variance et tests post-hoc")`

```{r assign_A10Ia_acidification_effect, echo=FALSE, results='asis'}
if (exists("assignment"))
  assignment("A10Ia_acidification_effect", part = NULL,
    url = "https://github.com/BioDataScience-Course/A10Ia_acidification_effect",
    course.ids = c(
      'S-BIOG-027' = !"A10Ia_{YY}M_acidification_effect"),
    course.urls = c(
      'S-BIOG-027' = "https://classroom.github.com/a/Jy9gGqx0"),
    course.starts = c(
      'S-BIOG-027' = !"{W[29]+1} 13:00:00"),
    course.ends = c(
      'S-BIOG-027' = !"{W[30]+1} 23:59:59"),
    term = "Q2", level = 3,
    toc = "ANOVA (effet de l'acidification)")
```

## Test de Kruskal-Wallis

Les conditions d'application de l'ANOVA sont assez restrictives, en particulier l'homoscédasticité et le normalité des résidus. Dans notre exemple, nous avons pu stabiliser la variance par une transformation puissance cinq, mais la distribution des résidus n'est pas optimale. Nous pouvons tout aussi bien décider de nous orienter vers la version non paramétrique équivalente : le **test de Kruskal-Wallis**.

```{block, type='warning'}
Le raisonnement entre ANOVA (test paramétrique) et Kruskal-Wallis (test non paramétrique) est le même qu'entre le test *t* de Student ou celui de Wilcoxon-Mann-Whitney. Lorsque les conditions sont remplies, l'ANOVA est un test plus puissant à utiliser en priorité. Nous le préférerons donc, sauf dans les cas impossibles où aucune transformation des données ne permet d'obtenir une distribution acceptable des résidus et une stabilisation des variances (homoscédasticité).
```

Le test de Kruskal-Wallis va comparer la localisation des points sur l'axe plutôt que les moyennes. Nous travaillons ici sur les **rangs**. La transformation en rangs consiste à ranger les observations de la plus petite à la plus grande et à remplacer les valeurs par leur position (leur rang) dans ce classement. Les ex-æquos reçoivent des rangs équivalents. Un petit exemple :

```{r}
# Un échantillon
x <- c(4.5, 2.1, 0.5, 2.4, 2.1, 3.5)
# Tri par ordre croissant
sort(x)
# Remplacement par les rangs
rank(sort(x))
```

Donc 0,5 donne 1, 2,1 apparaît deux fois en position 2 et 3, donc on leurs attribue à tous les deux le rang 2,5. 2,4 est remplacé par 4, 3,5 par 5 et enfin 4,5 est remplacé par 6.

Dans le test de Kruskal-Wallis, sous $H_0$ nous avons le même rang moyen (noté $mr$) pour chaque groupe parmi $k$. Sous $H_1$, au moins deux groupes ont des rangs moyens différents.

-   $H_0:mr_1 = mr_2 = ... = mr_k$

-   $H_1: \exists(i, j) \mathrm{\ tel\ que\ } mr_i \neq mr_j$

Nous ne détaillons pas les calculs (mais voyez [ici](http://vassarstats.net/textbook/ch14a.html)). Le calcul aboutit en fait à une valeur de $\chi^2_{obs}$ à $k - 1$ degrés de liberté lorsque *k* populations sont comparées. Ainsi, nous avons déjà un loi de distribution théorique pour le calcul de la valeur *P*. La zone de rejet est située à la droite de la distribution comme dans le cas de l'ANOVA et du test de $\chi^2$.

Le test est réalisé dans R avec la fonction `kruskal.test()`. Notez que les arguments sont les mêmes que pour l'ANOVA à un facteur.

<!-- Un "snippet" est disponible dans la SciViews Box (`.hn` pour `hypothesis tests: nonparametric`, ensuite `Kruskal-Wallis test`). -->

```{r}
kruskal.test(data = crabs2, aspect ~ group)
```

Ici, nous rejetons $H_0$ au seuil $\alpha$ de 5%. Nous pouvons dire qu'il y a au moins une différence significative entre les ratios au seuil $\alpha$ de 5% (test Kruskal-Wallis, $\chi^2$ = 140, ddl = 3, valeur *P* \<\< 10^-3^).

##### Conditions d'application {.unnumbered}

Le test de Kruskal-Wallis est naturellement moins contraignant que l'ANOVA.

-   échantillon représentatif (par exemple, aléatoire),
-   observations indépendantes,
-   variable dite **réponse** quantitative,
-   une variable dite **explicative** qualitative à trois niveaux ou plus,
-   les distributions au sein des différentes sous-populations sont, si possible, similaires, mais de forme quelconque.

### Test "post hoc" non paramétrique

Nous devrions également réaliser des tests "post hoc" lorsque nous rejetons $H_0$ du Kruskal-Wallis. Les versions non paramétriques des tests de comparaisons multiples sont moins connues. Nous pourrions effectuer des comparaisons deux à deux avec des tests de Wilcoxon-Mann-Whitney en appliquant une correction de Bonferroni.

Une autre option consiste à utiliser une version non paramétrique du test de Tukey basée sur les rangs, dit test post hoc non paramétrique asymptotique. Les explications assez techniques sont disponibles [ici](https://onlinelibrary.wiley.com/doi/abs/10.1002/1521-4036%28200109%2943%3A5%3C553%3A%3AAID-BIMJ553%3E3.0.CO%3B2-N). Ce test est disponible dans le package R {nparcomp}.

<!-- également disponible dans la SciViews Box depuis `.hn`, en sélectionnant l'entrée `Kruskal-Wallis: multiple comparisons` dans le menu déroulant. -->

```{r}
summary(crabs2_kw_comp <- nparcomp::nparcomp(data = crabs2, aspect ~ group))
plot(crabs2_kw_comp)
```

La présentation des résultats est plus détaillée que pour le HSD de Tukey. Pour la version textuelle, c'est la dernière colonne nommée `p.Value` qui nous intéresse et qui est à comparer au seuil $\alpha$ choisi. Le graphique est très similaire à la version paramétrique avec des barres d'erreurs autour de la différence des moyennes de chaque paire qui représente un intervalle de confiance à $1 - \alpha$. Ici, toutes les différences sont considérées comme significatives, même si la comparaison des mâles `B-M - O-M` avec une valeur *P* tout juste significative de 0,049 est à prendre avec des pincettes étant donné sa proximité du seuil. Notez que les intervalles de confiance ne sont plus symétriques dans le cas présent.

Au final, que conclure ? Lorsque l'ANOVA peut être utilisée, elle est à préférer. Ici avec une transformation puissance cinq, nous stabilisons la variance qui est le critère le plus sensible. De plus tant avec l'ANOVA qu'avec le Kruskal-Wallis, nous détectons des différences significatives. Les tests "post hoc" nous indiquent des différences significatives au seuil $\alpha$ de 5% sauf peut-être pour les mâles (il faudrait prévoir des mesures supplémentaires pour confirmer ou infirmer à ce niveau). Dans notre exemple, puisque les deux tests donnent sensiblement le même résultat, nous pouvons jouer la sécurité et présenter le Kruskal-Wallis dans notre publication. Ainsi, la question de la distribution des résidus ne se pose plus.

##### À vous de jouer ! {.unnumbered}

`r learnr("A10Lb_kruskal", title = "Test de Kruskal-Wallis", toc = "Test de Kruskal-Wallis")`

Appliquez vos nouvelles connaissances dans le projet de groupe que vous avez débuté aux deux modules précédents.

```{r assign_A10Ga_human_health_III, echo=FALSE, results='asis'}
if (exists("assignment2"))
  assignment2("A08Ga_human_health", part = "III",
    url = "https://github.com/BioDataScience-Course/A08Ga_human_health",
    course.ids = c(
      'S-BIOG-027' = !"A08Ga_{YY}M_human_health"),
    course.urls = c(
      'S-BIOG-027' = "https://classroom.github.com/a/HvBUqMqi"),
    course.starts = c(
      'S-BIOG-027' = !"{W[25]+1} 13:00:00"),
    course.ends = c(
      'S-BIOG-027' = !"{W[37]+1} 23:59:59"),
    term = "Q2", level = 4, n = 4,
    toc = "Biométrie humaine, part III")
```

```{=html}
<!-- **Note : Le projet suivant est la continuité d'un exercice débuté au premier quadrimestre et approfondis aux deux modules précédents.**

{r assign_A10Ga_urchin_IV, echo=FALSE, results='asis'}
if (exists("assignment2"))
  assignment2("A10Ga_urchin", part = "IV",
    url = "https://github.com/BioDataScience-Course/A03Ga_urchin",
    course.ids = c(
      'S-BIOG-027' = !"A03Ga_{YY}M_urchin", # TODO: how to account it in Q2???
      'S-BIOG-921' = !"A03Ga_{YY}C_urchin"),
    course.urls = c(
      'S-BIOG-027' = "https://classroom.github.com/a/...", 
      'S-BIOG-921' = "https://classroom.github.com/g/..."),
    course.starts = c(
      'S-BIOG-027' = !"{W[29]+1} 13:00:00"),
    course.ends = c(
      'S-BIOG-027' = !"{W[37]+1} 23:59:59"),
    term = "Q2", level = 4, n = 2,
    toc = "Biométrie de l'oursin violet, part IV") -->
```

## Récapitulatif des exercices

Vous venez de terminer le module 10. Ce module vous a permis d'apprendre une grande partie des subtilités qui se cachent derrière l'analyse de variance. Pour évaluer votre compréhension de cette matière, vous aviez les exercices suivants à réaliser :

`r show_ex_toc()`

##### Progression {.unnumbered}

`r launch_report("10", height = 800)`

````{=html}
<!--- G.E. Gardons nous cet exercice ?  
## Sciences des données et littérature

Vous avez à présent vu l'utilisation de plusieurs tests statistiques très courant en science (t-test, Anova, $\chi^2$,...). La restitution correcte de ce test au sein d'un publication scientifique est très importante.

```{block, type='bdd'}

Intéressez-vous à la restitution de résultat statistique au sein de la littérature. Un projet individuel est mis à votre disposition 

- <https://classroom.github.com/g/bc3Lh7C6>

```
--->
````