06-0_biostats.Rmd

---
site: bookdown::bookdown_site
output: bookdown::gitbook
---

# Biostatistiques {#chapitre-biostats}

 ***
️\ **Objectifs spécifiques**:

À la fin de ce chapitre, vous

- serez en mesure de définir les concepts de base en statistique: population, échantillon, variable, probabilité et distribution
- serez en mesure de calculer des statistiques descriptives de base: moyenne et écart-type, quartiles, maximum et minimum
- comprendrez les notions de test d'hypothèse, d'effet et de p-value, ainsi qu'éviter les erreurs communes dans leur interprétation
- saurez effectuer une modélisation statistique linéaire simple, multiple et mixte, entre autre sur des catégories
- saurez effectuer une modélisation statistique non linéaire simple, multiple et mixte

 ***

Aux chapitres précédents, nous avons vu comment visualiser, organiser et manipuler des tableaux de données. La statistique est une collection de disciplines liées à la collecte, l’organisation, l'analyse, l'interprétation et la présentation de données. Les biostatistiques sont des applications de ces disciplines à la biosphère.

Dans [*Principles and procedures of statistics: A biometrical approach*](https://www.amazon.com/Principles-Procedures-Statistics-Biometrical-Approach/dp/0070610282), Steel, Torie et Dickey (1997) définissent les statistiques ainsi:

> Les statistiques forment la science, pure et appliquée, de la création, du développement, et de l'application de techniques par lesquelles l'incertitude de l'induction inférentielle peut être évaluée. (ma traduction)

Alors que l'**inférence** consiste à généraliser des échantillons à l'ensemble d'une population, l'**induction** est un type de raisonnement qui permet de généraliser des observations sous forme de théories. En d'autres mots, les statistiques permettent d'évaluer l'incertitude sur des processus, de passer par infrence de l'échantillon à la population, puis par induction de passer de cette représentation d'une population en lois générales la concernant.

La définition de Whitlock et Schuluter (2015), dans [The Analysis of Biological Data](http://whitlockschluter.zoology.ubc.ca/), est plus simple et n'insiste que sur l'inférence:

> La statistique est l’étude des méthodes pour mesurer des aspects de populations à partir d’échantillons et pour quantifier l'incertitude des mesures. (ma traduction)

Les statistiques consistent à *faire du sens* (anglicisme assumé) avec des observations dans l'objectif de répondre à une question que vous aurez formulée clairement, préalablement à votre expérience.

<blockquote class="twitter-tweet" data-lang="fr"><p lang="en" dir="ltr">The more time I spend as The Statistician in the room, the more I think the best skill you can cultivate is the ability to remain calm and repeatedly ask &quot;What question are you trying to answer?&quot;</p>&mdash; Bryan Howie (@bryan_howie) <a href="https://twitter.com/bryan_howie/status/1073054519808876544?ref_src=twsrc%5Etfw">13 décembre 2018</a></blockquote>
<script async src="https://platform.twitter.com/widgets.js" charset="utf-8"></script>

Le flux de travail conventionnel en statistiques consiste à collecter des échantillons, transformer (prétraiter) les données, effectuer des tests, analyser les résultats, les interpréter et les visualiser.

$$Collecte \rightarrow Prétraitement \rightarrow Tests~statistiques \rightarrow Analyse \rightarrow Interprétation \rightarrow Visualisation $$

Ce chapitre à lui seul est trop court pour permettre d'intégrer toutes les connaissances nécessaires à une utilisation raisonnée des statistiques, mais fourni les bases pour aller plus loin. Notez que les erreurs d'interprétation statistiques sont courantes et la consultation de spécialistes n'est souvent pas un luxe. Mais bien que les statistiques soient complexes, *la plupart des opérations statistiques peuvent être effectuées sans l'assistance de statisticien.ne.s*... à condition de comprendre suffisamment les concepts utilisés.

Dans ce chapitre, nous verrons comment répondre correctement à une question valide et adéquate avec l'aide d'outils de calcul scientifique. Nous couvrirons les notions de bases des distributions et des variables aléatoires qui nous permettront d'effectuer des tests statistiques communs avec R. Nous couvrirons aussi les erreurs communément commises en recherche académique et les moyens simples de les éviter.

En plus des modules de base de R nous utiliserons

* les modules de la **`tidyverse`**,
* le module de données agricoles **`agridat`**, ainsi que
* le module **`nlme`** spécialisé pour la modélisation mixte.

Avant de survoler les applications statistiques avec R, je vais rapidement présenter quelques notions importantes en statistiques : populations et échantillons, variables, probabilités et distributions. Puis nous allons effectuer des tests d'hypothèse univariés (notamment les tests de *t* et les analyses de variance) et détailler la notion controversée de *p-value*. Je vais m'attarder plus longuement aux modèles linéaires généralisés, incluant en particulier des effets fixes et aléatoires (modèles mixtes), qui fournissent une trousse d'analyse polyvalente en analyse multivariée. Je terminerai avec les perspectives multivariées que sont les matrices de covariance et de corrélation.

## Populations et échantillons

Le principe d'inférence consiste à généraliser des conclusions à l'échelle d'une population à partir d'échantillons issus de cette population. Alors qu'une **population** contient tous les éléments étudiés, un **échantillon** d'une population est une observation unique. Une expérience bien conçue fera en sorte que les échantillons soient représentatifs de la population qui, la plupart du temps, ne peut être observée entièrement pour des raisons pratiques.

Les principes d'expérimentation servant de base à la conception d'une bonne méthodologie sont présentés dans le cours [*Dispositifs expérimentaux (BVG-7002)*](https://www.ulaval.ca/les-etudes/cours/repertoire/detailsCours/bvg-7002-dispositifs-experimentaux.html). Également, je recommande le livre *Principes d'expérimentation: planification des expériences et analyse de leurs résultats* de Pierre Dagnelie (2012), [disponible en ligne en format PDF](http://www.dagnelie.be/docpdf/ex2012.pdf). Un bon aperçu des dispositifs expérimentaux est aussi présenté dans [*Introductory Statistics with R*](https://www.springer.com/us/book/9780387790534), de Peter Dalgaard (2008), que vous pouvez télécharger [du site de la bibliothèque de l'Université Laval](https://ariane25.bibl.ulaval.ca/ariane/wicket/detail?c=ariane&m=S&rq.ct=PE&rq.fa=false&rq.r.esc=false&rq.r.l%5B0%5D.c=*&rq.r.l%5B0%5D.ex=false&rq.r.l%5B0%5D.op=AND&rq.r.l%5B0%5D.v=introductory+statistics+with+R&rq.r.la=*&rq.r.loc=*&rq.r.pft=false&rq.r.ta=*&rq.r.td=*&rq.rows=15&rq.st=0) vous avez un identifiant autorisé.

Une population est échantillonnée pour induire des **paramètres**: un rendement typique dans des conditions météorologiques, édaphiques et managériales données, la masse typique des faucons pèlerins, mâles et femelles, le microbiome typique d'un sol agricole ou forestier, etc. Une **statistique** est une estimation d'un paramètre calculée à partir des données, par exemple une moyenne et un écart-type, ou un intercept et une pente.

Par exemple, la moyenne ($\mu$) et l'écart-type ($\sigma$) d'une population sont estimés par les moyennes ($\bar{x}$) et écarts-types ($s$) calculés sur les données issues de l'échantillonnage.

Chaque paramètre est liée à une perspective que l'on désire connaître chez une population. Ces angles d'observations sont les **variables**.

## Les variables

Nous avons abordé au chapitre \@ref(chapitre-tableaux) la notion de *variable* par l'intermédiaire d'une donnée. Une variable est l'observation d'une caractéristique décrivant un échantillon. Si la charactéristique varie d'un échantillon à un autre sans que vous en expliquiez la raison (i.e. si identifier la source de la variabilité ne fait pas partie de votre expérience), on parlera de variable aléatoire. Même le hasard est régi par certaines lois: ce qui est aléatoire dans une variable peut être décrit par des **lois de probabilité**, que nous verrons plus bas.

Mais restons aux variables pour l'instant. Par convention, on peut attribuer aux variables un symbole mathématique. Par exemple, on peut donner à la masse volumique d'un sol (qui est le résultat d'une méthodologie précise) le symbole $\rho$. Lorsque l'on attribue une valeur à $\rho$, on parle d'une donnée. Chaque donnée d'une observation a un indice qui lui est propre, que l'on désigne souvent par $i$, que l'on place en indice $\rho_i$. Pour la première donnée, on a $i=1$, donc $\rho_1$. Pour un nombre $n$ d'échantillons, on aura $\rho_1$, $\rho_2$, $\rho_3$, ..., $\rho_n$, formant le vecteur $\rho = \left[\rho_1, \rho_2, \rho_3, ..., \rho_n \right]$.

En R, une variable est associée à un vecteur ou une colonne d'un tableau.

```{r biostats-data}
rho <- c(1.34, 1.52, 1.26, 1.43, 1.39) # matrice 1D
data <- data.frame(rho = rho) # tableau
data
```

Il existe plusieurs types de variables, qui se regroupe en deux grandes catégories: les **variables quantitatives** et les **variables qualitatives**.

### Variables quantitatives

Ces variables peuvent être continues dans un espace échantillonnal réel ou discrètes dans un espace échantillonnal ne considérant que des valeurs fixes. Notons que la notion de nombre réel est toujours une approximation en sciences expérimentales comme en calcul numérique, étant donnée que l'on est limité par la précision des appareils comme par le nombre d'octets à utiliser. Bien que les valeurs fixes des distributions discrètes ne soient pas toujours des valeurs entières, c'est bien souvent le cas en biostatistiques comme en démographie, où les décomptes d'individus sont souvent présents (et où la notion de fraction d'individus n'est pas acceptée).

### Variables qualitatives

On exprime parfois qu'une variable qualitative est une variable impossible à mesurer numériquement: une couleur, l'appartenance à espèce ou à une série de sol. Pourtant, dans bien des cas, les variables qualitatives peut être encodées en variables quantitatives. Par exemple, on peut accoler des pourcentages de sable, limon et argile à un loam sableux, qui autrement est décrit par la classe texturale d'un sol. Pour une couleur, on peut lui associer une longueur d'onde ou des pourcentages de rouge, vert et bleu, ainsi qu'un ton. En ce qui a trait aux variables ordonnées, il est possible de supposer un étalement. Par exemple, une variable d'intensité faible-moyenne-forte peut être transformée linéairement en valeurs quantitatives -1, 0 et 1. Attention toutefois, l'étalement peut parfois être quadratique ou logarithmique. Les séries de sol peuvent être encodées par la proportion de gleyfication ([Parent et al., 2017](https://www.frontiersin.org/articles/10.3389/fenvs.2017.00081/full#B4)). Quant aux catégories difficilement transformables en quantités, on pourra passer par l'**encodage catégoriel**, souvent appelé *dummyfication*, qui nous verrons plus loin. L'analyse qualitative consiste en l'analyse de verbatims, essentiellement utile en sciences sociales: nous n'en n'aurons pas besoin ici. Nous considérerons les variables qualitatives comme des variables quantitatives qui n'ont pas subi de prétraitement.

## Les probabilités

> « Nous sommes si éloignés de connaître tous les agens de la nature, et leurs divers modes d'action ; qu'il ne serait pas philosophique de nier les phénomènes, uniquement parce qu'ils sont inexplicables dans l'état actuel de nos connaissances. Seulement, nous devons les examiner avec une attention d'autant plus scrupuleuse, qu'il paraît plus difficile de les admettre ; et c'est ici que le calcul des probabilités devient indispensable, pour déterminer jusqu'à quel point il faut multiplier les observations ou les expériences, afin d'obtenir en faveur des agens qu'elles indiquent, une probabilité supérieure aux raisons que l'on peut avoir d'ailleurs, de ne pas les admettre. » — Pierre-Simon de Laplace

Une probabilité est la vraisemblance qu'un évènement se réalise chez un échantillon. Les probabilités forment le cadre des systèmes stochastiques, c'est-à-dire des systèmes trop complexes pour en connaître exactement les aboutissants, auxquels on attribue une part de hasard. Ces systèmes sont prédominants dans les processus vivants.

On peut dégager deux perspectives sur les probabilités: l'une passe par une interprétation fréquentielle, l'autre bayésienne.

- L'interprétation **fréquentielle** représente la fréquence des occurrences après un nombre infini d’évènements. Par exemple, si vous jouez à pile ou face un grand nombre de fois, le nombre de pile sera égal à la moitié du nombre de lancés. *L'approche fréquentielle teste si les données concordent avec un modèle du réel*. Il s'agit de l'interprétation communément utilisée.
- L'interprétation **bayésienne** vise à quantifier l'incertitude des phénomènes. Dans cette perspective, plus l'information s'accumule, plus l'incertitude diminue. Cette approche gagne en notoriété notamment parce qu'elle permet de décrire des phénomènes qui, intrinsèquement, ne peuvent être répétés infiniment (absence d'asymptote), comme celles qui sont bien définis dans le temps ou sur des populations limités. *L'approche bayésienne évalue la probabilité que le modèle soit réel*.

Une erreur courante consiste à aborder des statistiques fréquentielles comme des statistiques bayésiennes. Par exemple, si l'on désire évaluer la probabilité de l’existence de vie sur Mars, on devra passer par le bayésien, car avec les stats fréquentielles, l'on devra plutôt conclure si les données sont conformes ou non avec l'hypothèse de la vie sur Mars (exemple tirée du blogue [Dynamic Ecology](https://dynamicecology.wordpress.com/2011/10/11/frequentist-vs-bayesian-statistics-resources-to-help-you-choose/)).

Des rivalités factices s'installent enter les tenants des différentes approches, dont chacune, en réalité, répond à des questions différentes dont il convient réfléchir sur les limitations. Bien que les statistiques bayésiennes soient de plus en plus utilisées, nous ne couvrirons dans ce chapitre que l'approche fréquentielle. L'approche bayésienne est néanmoins traitée dans le chapitre \@ref(chapitre-biostats-bayes), qui est facultatif au cours.

## Les distributions

Une variable aléatoire peut prendre des valeurs selon des modèles de distribution des probabilités. Une distribution est une fonction mathématique décrivant la probabilité d'observer une série d'évènements. Ces évènements peuvent être des valeurs continues, des nombres entiers, des catégories, des valeurs booléennes (Vrai/Faux), etc. Dépendemment du type de valeur et des observations obtenues, on peut associer des variables à différentes lois de probabilité. Toujours, l'aire sous la courbe d'une distribution de probabilité est égale à 1.

En statistiques inférentielles, les distributions sont les modèles, comprenant certains paramètres comme la moyenne et la variance pour les distributions normales, à partir desquelles les données sont générées.

Il existe deux grandes familles de distribution: **discrètes** et **continues**. Les distributions discrètes sont contraintes à des valeurs prédéfinies (finies ou infinies), alors que les distributions continues prennent nécessairement un nombre infini de valeur, dont la probabilité ne peut pas être évaluée ponctuellement, mais sur un intervalle.

L'**espérance** mathématique est une fonction de tendance centrale, souvent décrite par un paramètre. Il s'agit de la moyenne d'une population pour une distribution normale. La **variance**, quant à elle, décrit la variabilité d'une population, i.e. son étalement autour de l'espérance. Pour une distribution normale, la variance d'une population est aussi appelée variance, souvent présentée par l'écart-type (égal à la racine carrée de la variance).

### Distribution binomiale

En tant que scénario à deux issues possibles, des tirages à pile ou face suivent une loi binomiale, comme toute variable booléenne prenant une valeur vraie ou fausse. En biostatistiques, les cas communs sont la présence/absence d'une espèce, d'une maladie, d'un trait phylogénétique, ainsi que les catégories encodées. Lorsque l'opération ne comprend qu'un seul échantillon (i.e. un seul tirage à pile ou face), il s'agit d'un cas particulier d'une loi binomiale que l'on nomme une loi de *Bernouilli*.

Pour 25 tirages à pile ou face indépendants (i.e. dont l'ordre des tirages ne compte pas), on peut dessiner une courbe de distribution dont la somme des probabilités est de 1. La fonction `dbinom` est une fonction de distribution de probabilités. Les fonctions de distribution de probabilités discrètes sont appelées des fonctions de masse.

```{r biostats-binom}
library("tidyverse")
x <- 0:25
y <- dbinom(x = x, size = 25, prob = 0.5)
print(paste('La somme des probabilités est de', sum(y)))
ggplot(data = tibble(x, y), mapping = aes(x, y)) +
  geom_segment(aes(x = x, xend = x, y = 0, yend = y), color = "grey50") +
  geom_point()
```

### Distribution de Poisson

La loi de Poisson (avec un P majuscule, introduite par le mathématicien français Siméon Denis Poisson et non pas l'animal) décrit des distributions discrètes de probabilité d'un nombre d’évènements se produisant dans l'espace ou dans le temps. Les distributions de Poisson décrivent ce qui tient du décompte. Il peut s'agir du nombre de grenouilles traversant une rue quotidiennement, du nombre de plants d'asclépiades se trouvant sur une terre cultivée, ou du nombre d’évènements de précipitation au mois de juin, etc. La distribution de Poisson n'a qu'un seul paramètre, $\lambda$, qui décrit tant la moyenne des décomptes.

Par exemple, en un mois de 30 jours, et une moyenne de 8 évènements de précipitation pour ce mois, on obtient la distribution suivante.

```{r biostats-poisson}
x <- 1:30
y <- dpois(x, lambda = 8)
print(paste('La somme des probabilités est de', sum(y)))
ggplot(data = data.frame(x, y), mapping = aes(x, y)) +
  geom_segment(aes(x = x, xend = x, y = 0, yend = y), color = "grey50") +
  geom_point()
```

### Distribution uniforme

La distribution la plus simple est probablement la distribution uniforme. Si la variable est discrète, chaque catégorie est associée à une probabilité égale. Si la variable est continue, la probabilité est directement proportionnelle à la largeur de l'intervalle. On utilise rarement la distribution uniforme en biostatistiques, sinon pour décrire des *a priori* vagues pour l'analyse bayésienne (ce sujet est traité dans le chapitre \@ref(chapitre-biostats-bayes)). Nous utilisons la fonction `dunif`. À la différence des distributions discrètes, les fonctions de distribution de probabilités continues sont appelées des fonctions de densité d'une loi de probabilité (*probability density function*).

```{r biostats-unif}
increment <- 0.01
x <- seq(-4, 4, by = increment)
y1 <- dunif(x, min = -3, max = 3)
y2 <- dunif(x, min = -2, max = 2)
y3 <- dunif(x, min = -1, max = 1)

print(paste('La somme des probabilités est de', sum(y3 * increment)))

gg_unif <- data.frame(x, y1, y2, y3) %>% gather(variable, value, -x)

ggplot(data = gg_unif, mapping = aes(x = x, y = value)) +
  geom_line(aes(colour = variable))
```

### Distribution normale

La plus répandue de ces lois est probablement la loi normale, parfois nommée loi gaussienne et plus rarement loi laplacienne. Il s'agit de la distribution classique en forme de cloche.

La loi normale est décrite par une moyenne, qui désigne la tendance centrale, et une variance, qui désigne l'étalement des probabilités autour de la moyenne. La racine carrée de la variance est l'écart-type.

Les distributions de mesures exclusivement positives (comme le poids ou la taille) sont parfois avantageusement approximées par une loi **log-normale**, qui est une loi normale sur le logarithme des valeurs: la moyenne d'une loi log-normale est la moyenne géométrique.

```{r biostats-norm}
increment <- 0.01
x <- seq(-10, 10, by = increment)
y1 <- dnorm(x, mean = 0, sd = 1)
y2 <- dnorm(x, mean = 0, sd = 2)
y3 <- dnorm(x, mean = 0, sd = 3)

print(paste('La somme des probabilités est de', sum(y3 * increment)))

gg_norm <- data.frame(x, y1, y2, y3) %>% gather(variable, value, -x)

ggplot(data = gg_norm, mapping = aes(x = x, y = value)) +
  geom_line(aes(colour = variable))
```

Quelle est la probabilité d'obtenir le nombre 0 chez une observation continue distribuée normalement dont la moyenne est 0 et l'écart-type est de 1? Réponse: 0. La loi normale étant une distribution continue, les probabilités non-nulles ne peuvent être calculés que sur des intervalles. Par exemple, la probabilité de retrouver une valeur dans l'intervalle entre -1 et 2 est calculée en soustrayant la probabilité cumulée à -1 de la probabilité cumulée à 2.

```{r biostats-norm-prob}
increment <- 0.01
x <- seq(-5, 5, by = increment)
y <- dnorm(x, mean = 0, sd = 1)

prob_between <- c(-1, 2)

gg_norm <- data.frame(x, y)
gg_auc <- gg_norm %>%
  filter(x > prob_between[1], x < prob_between[2]) %>%
  rbind(c(prob_between[2], 0)) %>%
  rbind(c(prob_between[1], 0))

ggplot(data.frame(x, y), aes(x, y)) +
  geom_polygon(data = gg_auc, fill = '#71ad50') + # #71ad50 est un code de couleur format hexadécimal
  geom_line()

prob_norm_between <- pnorm(q = prob_between[2], mean = 0, sd = 1) - pnorm(q = prob_between[1], mean = 0, sd = 1)
print(paste("La probabilité d'obtenir un nombre entre",
            prob_between[1], "et",
            prob_between[2], "est d'environ",
            round(prob_norm_between, 2) * 100, "%"))
```

La courbe normale peut être utile pour évaluer la distribution d'une population. Par exemple, on peut calculer les limites de région sur la courbe normale qui contient 95% des valeurs possibles en tranchant 2.5% de part et d'autre de la moyenne. Il s'agit ainsi de l'intervalle de confiance sur la déviation de la distribution.

```{r biostats-norm-ci}
increment <- 0.01
x <- seq(-5, 5, by = increment)
y <- dnorm(x, mean = 0, sd = 1)

alpha <- 0.05
prob_between <- c(qnorm(p = alpha/2, mean = 0, sd = 1),
                  qnorm(p = 1 - alpha/2, mean = 0, sd = 1))

gg_norm <- data.frame(x, y)
gg_auc <- gg_norm %>%
  filter(x > prob_between[1], x < prob_between[2]) %>%
  rbind(c(prob_between[2], 0)) %>%
  rbind(c(prob_between[1], 0))

ggplot(data = data.frame(x, y), mapping = aes(x, y)) +
  geom_polygon(data = gg_auc, fill = '#71ad50') + # #71ad50 est un code de couleur format hexadécimal
  geom_line() +
  geom_text(data = data.frame(x = prob_between,
                              y = c(0, 0),
                              labels = round(prob_between, 2)),
            mapping = aes(label = labels))
```

On pourrait aussi être intéressé à l'intervalle de confiance sur la moyenne. En effet, la moyenne suit aussi une distribution normale, dont la tendance centrale est la moyenne de la distribution, et dont l'écart-type est noté *erreur standard*. On calcule cette erreur en divisant la variance par le nombre d'observation, ou en divisant l'écart-type par la racine carrée du nombre d'observations. Ainsi, pour 10 échantillons:

```{rbiostats-norm-mean-ci}
increment <- 0.01
x <- seq(-5, 5, by = increment)
y <- dnorm(x, mean = 0, sd = 1)

alpha <- 0.05
prob_between <- c(qnorm(p = alpha/2, mean = 0, sd = 1) / sqrt(10),
                  qnorm(p = 1 - alpha/2, mean = 0, sd = 1) / sqrt(10))

gg_norm <- data.frame(x, y)
gg_auc <- gg_norm %>%
  filter(x > prob_between[1], x < prob_between[2]) %>%
  rbind(c(prob_between[2], 0)) %>%
  rbind(c(prob_between[1], 0))

ggplot(data = data.frame(x, y), mapping = aes(x, y)) +
  geom_polygon(data = gg_auc, fill = '#71ad50') + # #71ad50 est un code de couleur format hexadécimal
  geom_line() +
  geom_text(data = data.frame(x = prob_between,
                              y = c(0, 0),
                              labels = round(prob_between, 2)),
            mapping = aes(label = labels))
```

## Statistiques descriptives

On a vu comment générer des statistiques sommaires en R avec la fonction `summary()`. Reprenons les données d'iris.

```{r biostats-summ-oros}
data("iris")
summary(iris)
```

Pour précisément effectuer une moyenne et un écart-type sur un vecteur, passons par les fonctions `mean()` et `sd()`.

```{r biostats-stats-iris}
mean(iris$Sepal.Length)
sd(iris$Sepal.Length)
```

Pour effectuer un sommaire de tableau piloté par une fonction, nous passons par la gamme de fonctions `summarise()`, de dplyr. Dans ce cas, avec `group_by()`, nous fragmentons le tableau par espèce pour effectuer un sommaire sur toutes les variables.

```{r biostats-stats-iris-tv}
iris %>%
  group_by(Species) %>%
  summarise_all(mean)
```

Vous pourriez être intéressé par les quartiles à 25, 50 et 75%. Mais la fonction `summarise()` n'autorise que les fonctions dont la sortie est d'un seul objet, alors faisons sorte que l'objet soit une liste - lorsque l'on imbrique une fonction `funs`, le tableau à insérer dans la fonction est indiqué par un `.`.

```{r biostats-iris-tv-quantile}
iris %>%
  group_by(Species) %>%
  summarise_all(list(q25 = ~ quantile(., probs = 0.25),
                     q50 = ~ quantile(., probs = 0.50),
                     q75 = ~ quantile(., probs = 0.75)))
```

En mode programmation classique de R, on pourra générer les quartiles à la pièce.

```{r biostats-iris-classic-quantile}
quantile(iris$Sepal.Length[iris$Species == 'setosa'])
quantile(iris$Sepal.Length[iris$Species == 'versicolor'])
quantile(iris$Sepal.Length[iris$Species == 'virginica'])
```

La fonction `table()` permettra d'obtenir des décomptes par catégorie, ici par plages de longueurs de sépales. Pour obtenir les proportions du nombre total, il s'agit d'encapsuler le tableau croisé dans la fonction `prop.table()`.

```{r biostats-table}
tableau_croise <- table(iris$Species,
                        cut(iris$Sepal.Length, breaks = quantile(iris$Sepal.Length)))
tableau_croise
```

```{r biostats-prop-table}
prop.table(tableau_croise)
```

## Tests d'hypothèses à un et deux échantillons

Un test d'hypothèse permet de décider si une hypothèse est confirmée ou rejetée à un seuil de probabilité prédéterminé.

Cette section est inspirée du chapitre 5 de [Dalgaard, 2008](https://www.springer.com/us/book/9780387790534).

----

**Information: l'hypothèse nulle**. Les tests d'hypothèse évaluent des *effets* statistiques (qui ne sont pas nécessairement des effets de causalité). L'effet à évaluer peut être celui d'un traitement, d'indicateurs météorologiques (e.g. précipitations totales, degré-jour, etc.), de techniques de gestion des paysages, etc. Une recherche est menée pour évaluer l'hypothèse que l'on retrouve des différences entre des unités expérimentales. Par convention, l'**hypothèse nulle** (écrite $H_0$) est l'hypothèse qu'il n'y ait pas d'effet (c'est l'hypothèse de l'avocat du diable 😈) à l'échelle de la population (et non pas à l'échelle de l'échantillon). À l'inverse, l'**hypothèse alternative** (écrite $H_1$) est l'hypothèse qu'il y ait un effet à l'échelle de la population.

----

À titre d'exercice en stats, on débute souvent par en testant si deux vecteurs de valeurs continues proviennent de populations à moyennes différentes ou si un vecteur de valeurs a été généré à partir d'une population ayant une moyenne donner. Dans cette section, nous utiliserons la fonction `t.test()` pour les tests de t et la fonction `wilcox.test()` pour les tests de Wilcoxon (aussi appelé de Mann-Whitney).

### Test de t à un seul échantillon

Nous devons assumer, pour ce test, que l'échantillon est recueillit d'une population dont la distribution est normale, $\mathcal{N} \sim \left( \mu, \sigma^2 \right)$, et que chaque échantillon est indépendant l'un de l'autre. L'hypothèse nulle est souvent celle de l'avocat du diable, que la moyenne soit égale à une valeur donnée (donc la différence entre la moyenne de la population et une moyenne donnée est de zéro): ici, que $\mu = \bar{x}$. L'erreur standard sur la moyenne (ESM) de l'échantillon, $\bar{x}$ est calculée comme suit.

$$ESM = \frac{s}{\sqrt{n}}$$

où $s$ est l'écart-type de l'échantillon et $n$ est le nombre d'échantillons.

Pour tester l'intervalle de confiance de l'échantillon, on multiplie l'ESM par l'aire sous la courbe de densité couvrant une certaine proportion de part et d'autre de l'échantillon. Pour un niveau de confiance de 95%, on retranche 2.5% de part et d'autre.

```{r biostats-norm-error}
set.seed(33746)
x <- rnorm(20, 16, 4)

level <-  0.95
alpha <- 1-level

x_bar <- mean(x)
s <- sd(x)
n <- length(x)

error <- qnorm(1 - alpha/2) * s / sqrt(n)
error
```

intervalle de confiance est l'erreur de par et d'autre de la moyenne.

```{r biostats-norm-error-ci}
c(x_bar - error, x_bar + error)
```

Si la moyenne de la population est de 16, un nombre qui se situe dans l'intervalle de confiance on accepte l'hypothèse nulle au seuil 0.05. Si le nombre d'échantillon est réduit (généralement < 30), on passera plutôt par une distribution de t, avec $n-1$ degrés de liberté.

```{r biostats-student-error-ci}
error <- qt(1 - alpha/2, n-1) * s / sqrt(n)
c(x_bar - error, x_bar + error)
```

Plus simplement, on pourra utiliser la fonction `t.test()` en spécifiant la moyenne de la population. Nous avons généré 20 données avec une moyenne de 16 et un écart-type de 4. Nous savons donc que la vraie moyenne de l'échantillon est de 16. Mais disons que nous testons l'hypothèse que ces données sont tirées d'une population dont la moyenne est 18 (et implicitement que sont écart-type est de 4).

```{r biostats-ttest}
t.test(x, mu = 18)
```

La fonction retourne la valeur de t (*t-value*), le nombre de degrés de liberté ($n-1 = 19$), une description de l'hypothèse alternative (`alternative hypothesis: true mean is not equal to 18`), ainsi que l'intervalle de confiance au niveau de 95%. Le test contient aussi la *p-value*. Bien que la *p-value* soit largement utilisée en science

----

#### Information: la *p-value*

La *p-value*, ou valeur-p ou p-valeur, est utilisée pour trancher si, oui ou non, un résultat est **significatif**. En langage scientifique, le mot significatif ne devrait être utilisé *que* lorsque l'on réfère à un test d'hypothèse statistique. Vous retrouverez des *p-values* partout en stats. Les *p-values* indiquent la probabilité que les données ait été échantillonnées d'une population où un effet est observable selon le modèle statistique utilisé.

> La p-value est la probabilité que les données aient été générées pour obtenir un effet équivalent ou plus prononcé si l'hypothèse nulle est vraie.

Une *p-value* élevée indique que le modèle appliqué à vos données concorde avec la conclusion que l'hypothèse nulle est vraie, et inversement si la *p-value* est faible. Le seuil arbitraire utilisée en écologie et en agriculture, comme dans plusieurs domaines, est de 0.05. L'utilisation d'un seuil est toutefois contestée **à raison**. Une enquête menée dans la littérature scientifiques a révélé que 49% des 791 articles étudiés interprétaient un effet non significatif comme un effet nul ([Amrhein et al., 2019](https://www.nature.com/articles/d41586-019-00857-9)). En effet, une catégorisation de la p-value avec un seuil de significativité brouille le jugement sur l'importance des effets et de leur incertitude. Les six principes de l'[American Statistical Association](https://phys.org/news/2016-03-american-statistical-association-statement-significance.html) guident l'interprétation des *p-values*. [ma traduction]

1. Les *p-values* indique l'ampleur de l’incompatibilité des données avec le modèle statistique
2. Les *p-values* ne mesurent pas la probabilité que l'hypothèse étudiée soit vraie, ni la probabilité que les données ont été générées uniquement par la chance.
3. Les conclusions scientifiques et décisions d'affaire ou politiques ne devraient pas être basées sur l'atteinte d'une *p-value* à un seuil spécifique.
4. Une inférence appropriée demande un rapport complet et transparent.
5. Une *p-value*, ou une signification statistique, ne mesure pas l'ampleur d'un effet ou l'importance d'un résultat.
6. En tant que tel, une *p-value* n'offre pas une bonne mesure des évidences d'un modèle ou d'une hypothèse.

----

Dans le cas précédent, la *p-value* était de 0.01014. Pour aider notre interprétation, prenons l'hypothèse alternative: `true mean is not equal to 18`. L'hypothèse nulle était bien que *la vraie moyenne est égale à 18*. Insérons la *p-value* dans la définition: la probabilité que les données aient été générées pour obtenir un effet équivalent ou plus prononcé si l'hypothèse nulle est vraie est de 0.01014. Il est donc très peu probable que les données soient tirées d'un échantillon dont la moyenne est de 18. Au seuil de signification de 0.05, on rejette l'hypothèse nulle et l'on conclut qu'à ce seuil de confiance, l'échantillon ne provient pas d'une population ayant une moyenne de 18.

### Attention: mauvaises interprétations des *p-values*

> "La p-value n'a jamais été conçue comme substitut au raisonnement scientifique" [Ron Wasserstein, directeur de l'American Statistical Association](https://phys.org/news/2016-03-american-statistical-association-statement-significance.html) [ma traduction].

**Un résultat montrant une p-value plus élevée que 0.05 est-il pertinent?**

Lors d'une conférence, Dr Evil ne présentent que les résultats significatifs de ses essais au seuil de 0.05. Certains essais ne sont pas significatifs, mais bon, ceux-ci ne sont pas importants... En écartant ces résultats, Dr Evil commet 3 erreurs:

1. La *p-value* n'est pas un bon indicateur de l'importance d'un test statistique. L'importance d'une variable dans un modèle devrait être évaluée par la valeur de son coefficient. Son incertitude devrait être évaluée par sa variance. Une manière d'évaluer plus intuitive la variance est l'écart-type ou l'intervalle de confiance. À un certain seuil d'intervalle de confiance, la p-value traduira la probabilité qu'un coefficient soit réellement nul ait pu générer des données démontrant un coefficient égal ou supérieur.
1. Il est tout aussi important de savoir que le traitement fonctionne que de savoir qu'il ne fonctionne pas. Les résultats démontrant des effets sont malheureusement davantage soumis aux journaux et davantage publiés que ceux ne démontrant pas d'effets ([Decullier et al., 2005]( https://doi.org/10.1136/bmj.38488.385995.8F )).
1. Le seuil de 0.05 est arbitraire.

----

#### Attention au *p-hacking*

Le *p-hacking* (ou *data dredging*) consiste à manipuler les données et les modèles pour faire en sorte d'obtenir des *p-values* favorables à l'hypothèse testée et, éventuellement, aux conclusions recherchées. **À éviter dans tous les cas. Toujours. Toujours. Toujours.**

```{r biostats-johnoliver, fig.align="center", fig.cap="Un sketch humoristique de John Oliver sur le *p-hacking*, [Last week tonight, 2016](https://youtu.be/0Rnq1NpHdmw) (en anglais)", echo = FALSE}
knitr::include_graphics("images/05_p-hacking.png")
```

----

### Test de Wilcoxon à un seul échantillon

Le test de t suppose que la distribution des données est normale... ce qui est rarement le cas, surtout lorsque les échantillons sont peu nombreux. Le test de Wilcoxon ne demande aucune supposition sur la distribution: c'est un test non-paramétrique basé sur le tri des valeurs.

```{r biostats-wilcox-test}
wilcox.test(x, mu = 18)
```

Le `V` est la somme des rangs positifs. Dans ce cas, la *p-value* est semblable à celle du test de t, et les mêmes conclusions s'appliquent.

### Tests de t à deux échantillons

Les tests à un échantillon servent plutôt à s'exercer: rarement en aura-t-on besoin en recherche, où plus souvent, on voudra comparer les moyennes de deux unités expérimentales. L'expérience comprend donc deux séries de données continues, $x_1$ et $x_2$, issus de lois de distribution normale $\mathcal{N} \left( \mu_1, \sigma_1^2 \right)$ et $\mathcal{N} \left( \mu_2, \sigma_2^2 \right)$, et nous testons l'hypothèse nulle que $\mu_1 = \mu_2$. La statistique t est calculée comme suit.

$$t = \frac{\bar{x_1} - \bar{x_2}}{ESDM}$$

L'ESDM est l'erreur standard de la différence des moyennes:

$$ESDM = \sqrt{ESM_1^2 + ESM_2^2}$$

Si vous supposez que les variances sont identiques, l'erreur standard (s) est calculée pour les échantillons des deux groupes, puis insérée dans le calcul des ESM. La statistique t sera alors évaluée à $n_1 + n_2 - 2$ degrés de liberté. Si vous supposez que la variance est différente (*procédure de Welch*), vous calculez les ESM avec les erreurs standards respectives, et la statistique t devient une approximation de la distribution de t avec un nombre de degrés de liberté calculé à partir des erreurs standards et du nombre d'échantillon dans les groupes: cette procédure est considérée comme plus prudente ([Dalgaard, 2008](https://www.springer.com/us/book/9780387790534), page 101).

Prenons les données d'iris pour l'exemple en excluant l'iris setosa étant donnée que les tests de t se restreignent à deux groupes. Nous allons tester la longueur des pétales.

```{r biostats-iris-petal}
iris_pl <- iris %>%
    filter(Species != "setosa") %>%
    select(Species, Petal.Length)
sample_n(iris_pl, 5)
```

Dans la prochaine cellule, nous introduisons l'*interface-formule* de R, où l'on retrouve typiquement le `~`, entre les variables de sortie à gauche et les variables d'entrée à droite. Dans notre cas, la variable de sortie est la variable testée, `Petal.Length`, qui varie en fonction du groupe `Species`, qui est la variable d'entrée (variable explicative) - nous verrons les types de variables plus en détails dans la section [Les modèles statistiques](#Les-mod%C3%A8les-statistiques), plus bas.

```{r biostats-iris-petal-ttest}
t.test(formula = Petal.Length ~ Species,
       data = iris_pl, var.equal = FALSE)
```

Nous obtenons une sortie similaire aux précédentes. L'intervalle de confiance à 95% exclu le zéro, ce qui est cohérent avec la p-value très faible, qui nous indique le rejet de l'hypothèse nulle au seuil 0.05. Les données montrent que les groupes ont des moyennes de longueurs de pétale différentes.

----

#### Enregistrer les résultats d'un test

Il est possible d'enregistrer un test dans un objet.

```{r biostats-save-test}
tt_pl <- t.test(formula = Petal.Length ~ Species,
                data = iris_pl, var.equal = FALSE)
summary(tt_pl)
str(tt_pl)
```

----

### Comparaison des variances

Pour comparer les variances, on a recours au test de F (F pour Fisher).

```{r biostats-ftest}
var.test(formula = Petal.Length ~ Species,
         data = iris_pl)
```

Il semble que l'on pourrait relancer le test de *t* sans la procédure Welch, avec `var.equal = TRUE`.

### Tests de Wilcoxon à deux échantillons

Cela ressemble au test de t!

```{r biostats-wilcox-2ech}
wilcox.test(formula = Petal.Length ~ Species,
       data = iris_pl, var.equal = TRUE)
```

### Les tests pairés

Les tests pairés sont utilisés lorsque deux échantillons proviennent d'une même unité expérimentale: il s'agit en fait de tests sur la différence entre deux observations.

```{r biostats-ap}
set.seed(2555)

n <- 20
avant <- rnorm(n, 16, 4)
apres <- rnorm(n, 18, 3)
```

Il est important de spécifier que le test est pairé, la valeur par défaut de `paired` étant `FALSE`.

```{r biostats-ttest-paired}
t.test(avant, apres, paired = TRUE)
```

L'hypothèse nulle qu'il n'y ait pas de différence entre l'avant et l'après traitement est acceptée au seuil 0.05.

**Exercice**. Effectuer un test de Wilcoxon pairé.

## L'analyse de variance

L'analyse de variance consiste à comparer des moyennes de plusieurs groupe distribués normalement et de même variance. Cette section sera élaborée prochainement plus en profondeur. Considérons-la pour le moment comme une régression sur une variable catégorielle.

```{r biostats-anova}
pl_aov <- aov(Petal.Length ~ Species, iris)
summary(pl_aov)
```

La prochaine section, justement, est vouée aux modèles statistiques explicatifs, qui incluent la régression.

## Les modèles statistiques

La modélisation statistique consiste à lier de manière explicite des variables de sortie $y$ (ou variables-réponse ou variables dépendantes) à des variables explicatives $x$ (ou variables prédictives / indépendantes / covariables). Les variables-réponse sont modélisées par une fonction des variables explicatives ou prédictives.

Pourquoi garder les termes *explicatives* et *prédictives*? Parce que les modèles statistiques (basés sur des données et non pas sur des mécanismes) sont de deux ordres. D'abord, les modèles **prédictifs** sont conçus pour prédire de manière fiable une ou plusieurs variables-réponse à partir des informations contenues dans les variables qui sont, dans ce cas, prédictives. Ces modèles sont couverts dans le chapitre 11 de ce manuel (en développement). Lorsque l'on désire tester des hypothèses pour évaluer quelles variables expliquent la réponse, on parlera de modélisation (et de variables) **explicatives**. En inférence statistique, on évaluera les *corrélations* entre les variables explicatives et les variables-réponse. Un lien de corrélation n'est pas un lien de causalité. L'inférence causale peut en revanche être évaluée par des [*modèles d'équations structurelles*](https://www.amazon.com/Cause-Correlation-Biology-Structural-Equations/dp/1107442591), sujet qui fera éventuellement partie de ce cours.

Cette section couvre la modélisation explicative. Les variables qui contribuent à créer les modèles peuvent être de différentes natures et distribuées selon différentes lois de probabilité. Alors que les modèles linéaires simples (*lm*) impliquent une variable-réponse distribuée de manière continue, les modèles linéaires généralisés peuvent aussi expliquer des variables de sorties discrètes.

Dans les deux cas, on distinguera les variables fixes et les variables aléatoires. Les **variables fixes** sont les variables testées lors de l'expérience: dose du traitement, espèce/cultivar, météo, etc. Les **variables aléatoires** sont les sources de variation qui génèrent du bruit dans le modèle: les unités expérimentales ou le temps lors de mesures répétées. Les modèles incluant des effets fixes seulement sont des modèles à effets fixes. Généralement, les modèles incluant des variables aléatoires incluent aussi des variables fixes: on parlera alors de modèles mixtes. Nous couvrirons ces deux types de modèle.

### Modèles à effets fixes

Les tests de t et de Wilcoxon, explorés précédemment, sont des modèles statistiques à une seule variable. Nous avons vu dans l'*interface-formule* qu'une variable-réponse peut être liée à une variable explicative avec le tilde `~`. En particulier, le test de t est régression linéaire univariée (à une seule variable explicative) dont la variable explicative comprend deux catégories. De même, l'anova est une régression linéaire univariée dont la variable explicative comprend plusieurs catégories. Or l'interface-formule peut être utilisé dans plusieurs circonstances, notamment pour ajouter plusieurs variables de différents types: on parlera de régression multivariée.

La plupart des modèles statistiques peuvent être approximés comme une combinaison linéaire de variables: ce sont des modèles linéaires. Les modèles non-linéaires impliquent des stratégies computationnelles complexes qui rendent leur utilisation plus difficile à manœuvrer.

Un modèle linéaire univarié prendra la forme $y = \beta_0 + \beta_1 x + \epsilon$, où $\beta_0$ est l'intercept et $\beta_1$ est la pente et $\epsilon$ est l'erreur.

Vous verrez parfois la notation $\hat{y} = \beta_0 + \beta_1 x$. La notation avec le chapeau $\hat{y}$ exprime qu'il s'agit des valeurs générées par le modèle. En fait, $y = \hat{y} - \epsilon$.

#### Modèle linéaire univarié avec variable continue

Prenons les données [`lasrosas.corn`](https://rdrr.io/cran/agridat/man/lasrosas.corn.html) incluses dans le module `agridat`, où l'on retrouve le rendement d'une production de maïs à dose d'azote variable, en Argentine.

```{r biostats-load-agridat}
library("agridat")
data("lasrosas.corn")
sample_n(lasrosas.corn, 10)
```

Ces données comprennent plusieurs variables. Prenons le rendement (`yield`) comme variable de sortie et, pour le moment, ne retenons que la dose d'azote (`nitro`) comme variable explicative: il s'agit d'une régression univariée. Les deux variables sont continues. Explorons d'abord le nuage de points de l'une et l'autre.

```{r biostats-plot-nitro-agridat}
ggplot(data = lasrosas.corn, mapping = aes(x = nitro, y = yield)) +
    geom_point()
```

L'hypothèse nulle est que la dose d'azote n'affecte pas le rendement, c'est à dire que le coefficient de pente et nul. Une autre hypothèse est que l'intercept est nul: donc qu'à dose de 0, rendement de 0. Un modèle linéaire à variable de sortie continue est créé avec la fonction `lm()`, pour *linear model*.

```{r biostats-lm1-agridat}
modlin_1 <- lm(yield ~ nitro, data = lasrosas.corn)
summary(modlin_1)
```

Le diagnostic du modèle comprend plusieurs informations. D'abord la formule utilisée, affichée pour la traçabilité. Viens ensuite un aperçu de la distribution des résidus. La médiane devrait s'approcher de la moyenne des résidus (qui est toujours de 0). Bien que le -3.079 peut sembler important, il faut prendre en considération de l'échelle de y, et ce -3.079 est exprimé en terme de rendement, ici en quintaux (i.e. 100 kg) par hectare. La distribution des résidus mérite d'être davantage investiguée. Nous verrons cela un peu plus tard.

Les coefficients apparaissent ensuite. Les estimés sont les valeurs des effets. R fournit aussi l'erreur standard associée, la valeur de t ainsi que la p-value (la probabilité d'obtenir cet effet ou un effet plus extrême si en réalité il y avait absence d'effet). L'intercept est bien sûr plus élevé que 0 (à dose nulle, on obtient 65.8 quintaux par hectare en moyenne). La pente de la variable `nitro` est de ~0.06: pour chaque augmentation d'un kg/ha de dose, on a obtenu ~0.06 quintaux/ha de plus de maïs. Donc pour 100 kg/ha de N, on a obtenu un rendement moyen de 6 quintaux de plus que l'intercept. Soulignons que l'ampleur du coefficient est très important pour guider la fertilisation: ne rapporter que la p-value, ou ne rapporter que le fait qu'elle est inférieure à 0.05 (ce qui arrive souvent dans la littérature), serait très insuffisant pour l'interprétation des statistiques. La p-value nous indique néanmoins qu'il serait très improbable qu'une telle pente ait été générée alors que celle-ci est nulle en réalité. Les étoiles à côté des p-values indiquent l'ampleur selon l'échelle `Signif. codes` indiquée en-dessous du tableau des coefficients.

Sous ce tableau, R offre d'autres statistiques. En outre, les R² et R² ajustés indiquent si la régression passe effectivement par les points. Le R² prend un maximum de 1 lorsque la droite passe exactement sur les points.

Enfin, le test de F génère une p-value indiquant la probabilité que les coefficients de pente ait été générés si les vrais coefficients étaient nuls. Dans le cas d'une régression univariée, cela répète l'information sur l'unique coefficient.

On pourra également obtenir les intervalles de confiance avec la fonction `confint()`.

```{r biostats-lm1-agridat-ci}
confint(modlin_1, level = 0.95)
```

Ou soutirer l'information de différentes manières, comme avec la fonction `coefficients()`.

```{r biostats-lm1-agridat-coeff}
coefficients(modlin_1)
```

Également, on pourra exécuter le modèle sur les données qui ont servi à le générer:

```{r biostats-lm1-agridat-predict}
predict(modlin_1)[1:5]
```

Ou sur des données externes.

```{r biostats-lm1-agridat-predict-newdata}
nouvelles_donnees <- data.frame(nitro = seq(from = 0, to = 100, by = 5))
predict(modlin_1, newdata = nouvelles_donnees)[1:5]
```

#### Analyse des résidus

Les résidus sont les erreurs du modèle. C'est le vecteur $\epsilon$, qui est un décalage entre les données et le modèle. Le R² est un indicateur de l'ampleur du décalage, mais une régression linéaire explicative en bonne et due forme devrait être accompagnée d'une analyse des résidus. On peut les calculer par $\epsilon = y - \hat{y}$, mais aussi bien utiliser la fonction `residuals()`.

```{r biostats-lm1-residual}
res_df <- data.frame(nitro = lasrosas.corn$nitro,
                     residus_lm = residuals(modlin_1),
                     residus_calcul = lasrosas.corn$yield - predict(modlin_1))
sample_n(res_df, 10)
```

Dans une bonne régression linéaire, on ne retrouvera pas de structure identifiable dans les résidus, c'est-à-dire que les résidus sont bien distribués de part et d'autre du modèle de régression.

```{r biostats-lm1-residuals-plot}
ggplot(res_df, aes(x = nitro, y = residus_lm)) +
  geom_point() +
  labs(x = "Dose N", y = "Résidus") +
  geom_hline(yintercept = 0, col = "red", size = 1)
```

Bien que le jugement soit subjectif, on peut dire avec confiance qu'il n'y a pas structure particulière. En revanche, on pourrait générer un $y$ qui varie de manière quadratique avec $x$, un modèle linéaire montrera une structure évidente.

```{r biostats-residuals-example-banana}
set.seed(36164)
x <- 0:100
y <- 10 + x*1 + x^2 * 0.05 + rnorm(length(x), 0, 50)
modlin_2 <- lm(y ~ x)
ggplot(data.frame(y, residus = residuals(modlin_2)),
       aes(x = x, y = residus)) +
  geom_point() +
  labs(x = "x", y = "Résidus") +
  geom_hline(yintercept = 0, col = "red", size = 1)
```
De même, les résidus ne devraient pas croître avec $x$.

```{r biostats-residuals-example-expand}
set.seed(3984)
x <- 0:100
y <-  10 + x + x * rnorm(length(x), 0, 2)
modlin_3 <- lm(y ~ x)
ggplot(data.frame(x, residus = residuals(modlin_3)),
       aes(x = x, y = residus)) +
  geom_point() +
  labs(x = "x", y = "Résidus") +
  geom_hline(yintercept = 0, col = "red", size = 1)
```

On pourra aussi inspecter les résidus avec un graphique de leur distribution. Reprenons notre modèle de rendement du maïs.

```{r biostats-residuals-agridat-groups}
ggplot(res_df, aes(x = residus_lm)) +
  geom_histogram(binwidth = 2, color = "white") +
  labs(x = "Residual")
```

L'histogramme devrait présenter une distribution normale. Les tests de normalité comme le test de Shapiro-Wilk peuvent aider, mais ils sont généralement très sévères.

```{r biostats-residuals-shapiro}
shapiro.test(res_df$residus_lm)
```

L'hypothèse nulle que la distribution est normale est rejetée au seuil 0.05. Dans notre cas, il est évident que la sévérité du test n'est pas en cause, car les résidus semble générer trois ensembles. Ceci indique que les variables explicatives sont insuffisantes pour expliquer la variabilité de la variable-réponse.

#### Régression multiple

Comme c'est le cas pour bien des phénomènes en écologie, le rendement d'une culture n'est certainement pas expliqué seulement par la dose d'azote.

Lorsque l'on combine plusieurs variables explicatives, on crée un modèle de régression multivariée, ou une régression multiple. Bien que les tendances puissent sembler non-linéaires, l'ajout de variables et le calcul des coefficients associés reste un problème d'algèbre linéaire.

On pourra en effet généraliser les modèles linéaires, univariés et multivariés, de la manière suivante.

$$ y = X \beta + \epsilon $$

où:

$X$ est la matrice du modèle à $n$ observations et $p$ variables.

$$ X = \left( \begin{matrix}
1 & x_{11} & \cdots & x_{1p}  \\
1 & x_{21} & \cdots & x_{2p}  \\
\vdots & \vdots & \ddots & \vdots  \\
1 & x_{n1} & \cdots & x_{np}
\end{matrix} \right) $$

$\beta$ est la matrice des $p$ coefficients, $\beta_0$ étant l'intercept qui multiplie la première colonne de la matrice $X$.

$$ \beta = \left( \begin{matrix}
\beta_0  \\
\beta_1  \\
\vdots \\
\beta_p
\end{matrix} \right) $$

$\epsilon$ est l'erreur de chaque observation.

$$ \epsilon = \left( \begin{matrix}
\epsilon_0  \\
\epsilon_1  \\
\vdots \\
\epsilon_n
\end{matrix} \right) $$

#### Modèles linéaires univariés avec variable catégorielle **nominale**

Une variable catégorielle nominale (non ordonnée) utilisée à elle seule dans un modèle comme variable explicative, est un cas particulier de régression multiple. En effet, l'**encodage catégoriel** (ou *dummyfication*) transforme une variable catégorielle nominale en une matrice de modèle comprenant une colonne désignant l'intercept (une série de 1) désignant la catégorie de référence, ainsi que des colonnes pour chacune des autres catégories désignant l'appartenance (1) ou la non appartenance (0) de la catégorie désignée par la colonne.

##### L'encodage catégoriel

Une variable à $C$ catégories pourra être déclinée en $C$ variables dont chaque colonne désigne par un 1 l'appartenance au groupe de la colonne et par un 0 la non-appartenance. Pour l'exemple, créons un vecteur désignant le cultivar de pomme de terre.

```{r biostats-dummy-coding}
data <- data.frame(cultivar = factor(c('Superior', 'Superior', 'Superior', 'Russet', 'Kenebec', 'Russet')))
model.matrix(~cultivar, data)
```

Nous avons trois catégories, encodées en trois colonnes. La première colonne est un intercept et les deux autres décrivent l'absence (0) ou la présence (1) des cultivars Russet et Superior. Le cultivar Kenebec est absent du tableau. En effet, en partant du principe que l'appartenance à une catégorie est mutuellement exclusive, c'est-à-dire qu'un échantillon ne peut être assigné qu'à une seule catégorie, on peut déduire une catégorie à partir de l'information sur toutes les autres. Par exemple, si `cultivar_Russet` et `cultivar_Superior` sont toutes deux égales à $0$, on conclura que `cultivar_Kenebec` est nécessairement égal à $1$. Et si l'un d'entre `cultivar_Russet` et `cultivar_Superior` est égal à $1$, `cultivar_Kenebec` est nécessairement égal à $0$. L'information contenue dans un nombre $C$ de catégorie peut être encodée dans un nombre $C-1$ de colonnes. C'est pourquoi, dans une analyse statistique, on désignera une catégorie comme une référence, que l'on détecte lorsque toutes les autres catégories sont encodées avec des $0$: cette référence sera incluse dans l'intercept. La catégorie de référence par défaut en R est celle la première catégorie dans l'ordre alphabétique. On pourra modifier cette référence avec la fonction `relevel()`.

```{r biostats-dummy-coding-relevel}
data$cultivar <- relevel(data$cultivar, ref = "Superior")
model.matrix(~cultivar, data)
```

Pour certains modèles, vous devrez vous assurer vous-même de l'encodage catégoriel. Pour d'autre, en particulier avec l'*interface par formule* de R, ce sera fait automatiquement.

##### Exemple d'application

Prenons la topographie du terrain, qui peut prendre plusieurs niveaux.

```{r biostats-lasrosas-topo-levels}
levels(lasrosas.corn$topo)
```

Explorons le rendement selon la topographie.

```{r biostats-lasrosas-topo-boxplot}
ggplot(lasrosas.corn, aes(x = topo, y = yield)) +
    geom_boxplot()
```

Les différences sont évidentes, et la modélisation devrait montrer des effets différents.

L'encodage catégoriel peut être visualisé en générant la matrice de modèle avec la fonction `model.matrix()` et l'interface-formule - sans la variable-réponse.

```{r biostats-lasrosas-topo-model-matrix}
model.matrix(~ topo, data = lasrosas.corn) %>%
    tbl_df() %>% # tbl_df pour transformer la matrice en tableau
    sample_n(10)
```

Dans le cas d'un modèle avec une variable catégorielle nominale seule, l'intercept représente la catégorie de référence, ici `E`. Les autres colonnes spécifient l'appartenance (1) ou la non-appartenance (0) de la catégorie pour chaque observation.

Cette matrice de modèle utilisée pour la régression donnera un intercept, qui indiquera l'effet de la catégorie de référence, puis les différences entre les catégories subséquentes et la catégorie de référence.

```{r biostats-lasrosas-topo-levels-lm}
modlin_4 <- lm(yield ~ topo, data = lasrosas.corn)
summary(modlin_4)
```

Le modèle linéaire est équivalent à l'anova, mais les résultats de `lm` sont plus élaborés.

```{r biostats-lasrosas-topo-levels-summ}
summary(aov(yield ~ topo, data = lasrosas.corn))
```

L'analyse de résidus peut être effectuée de la même manière.

#### Modèles linéaires univariés avec variable catégorielle **ordinale**

Bien que j'introduise la régression sur variable catégorielle ordinale à la suite de la section sur les variables nominales, nous revenons dans ce cas à une régression simple, univariée. Voyons un cas à 5 niveaux.

```{r biostats-encoding-ordinal}
statut <- c("Totalement en désaccord",
            "En désaccord",
            "Ni en accord, ni en désaccord",
            "En accord",
            "Totalement en accord")
statut_o <- factor(statut, levels = statut, ordered=TRUE)
model.matrix(~statut_o) # ou bien, sans passer par model.matrix, contr.poly(5) où 5 est le nombre de niveaux
```

La matrice de modèle a 5 colonnes, soit le nombre de niveaux: un intercept, puis 4 autres désignant différentes valeurs que peuvent prendre les niveaux. Ces niveaux croient-ils linéairement? De manière quadratique, cubique ou plus loin dans des distributions polynomiales?

```{r biostats-encoding-ordinal-plot}
modmat_tidy <- data.frame(statut, model.matrix(~statut_o)[, -1]) %>%
    gather(variable, valeur, -statut)
modmat_tidy$statut <- factor(modmat_tidy$statut,
                             levels = statut,
                             ordered=TRUE)
ggplot(data = modmat_tidy, mapping = aes(x = statut, y = valeur)) +
    facet_wrap(. ~ variable) +
    geom_point() +
    theme(axis.text.x = element_text(angle = 45, hjust = 1))
```

Règle générale, pour les variables ordinales, on préférera une distribution linéaire, et c'est l'option par défaut de la fonction `lm()`. L'utilisation d'une autre distribution peut être effectuée à la mitaine en utilisant dans le modèle la colonne désirée de la sortie de la fonction `model.matrix()`.

#### Régression multiple à plusieurs variables

Reprenons le tableau de données du rendement de maïs.

```{r biostats-multivariate-lm-data}
head(lasrosas.corn)
```

Pour ajouter des variables au modèle dans l'interface-formule, on additionne les noms de colonne. La variable `lat` désigne la latitude, la variable `long` désigne la latitude et la variable `bv` (*brightness value*) désigne la teneur en matière organique du sol (plus `bv` est élevée, plus faible est la teneur en matière organique).

```{r biostats-multivariate-lm}
modlin_5 <- lm(yield ~ lat + long + nitro + topo + bv,
               data = lasrosas.corn)
summary(modlin_5)
```

L'ampleur des coefficients est relatif à l'échelle de la variable. En effet, un coefficient de 5541 sur la variable `lat` n'est pas comparable au coefficient de la variable `bv`, de -0.5089, étant donné que les variables ne sont pas exprimées avec la même échelle. Pour les comparer sur une même base, on peut centrer (soustraire la moyenne) et réduire (diviser par l'écart-type).

```{r biostats-multivariate-lm-scaled}
lasrosas.corn_sc <- lasrosas.corn %>%
    mutate_at(c("lat", "long", "nitro", "bv"), scale)

modlin_5_sc <- lm(yield ~ lat + long + nitro + topo + bv,
               data = lasrosas.corn_sc)
summary(modlin_5_sc)
```

Typiquement, les variables catégorielles, qui ne sont pas mises à l'échelle, donneront des coefficients plus élevées, et devrons être évaluées entre elles et non comparativement aux variables mises à l'échelle. Une manière conviviale de représenter des coefficients consiste à utiliser la fonction `tidy` du module **`broom`**, qui génère un tableau contennt les coefficients ainsi que leurs intervalles de confiance, que nous pourrons ensuite porter graphiquement.

```{r biostats-multivariate-lm-ci}
library("broom") # ou bien charger le méta-module tidymodels
intervals <- tidy(modlin_5_sc, conf.int = TRUE, conf.level = 0.95)
intervals
```

La valeur par défaut de l'argument `conf.level` est de 0.95, mais je vous suggère de toujours l'écrire de manière explicite, ne serait-ce que pour rappeler à vous-même ainsi qu'à vos collègues, que cette valeur est arbitraire: il s'agit d'une décision d'analyse, non pas d'une valeur à utiliser par convention.

Pour le graphique, on aura avantage à séparer les effets catégoriels aux effets numériques pour mieux interpréter leurs effets entre eux. J'utilise la fonction `dplyr::case_when()` pour créer une nouvelle colonne qui catégorisera les termes de l'équation. Cette catégorie me permettra d'effectuer un `facet_wrap()`.

```{r biostats-multivariate-lm-ci-plot}
intervals %>% 
  mutate(type = case_when(
    term %in% c("topoHT", "topoLO", "topoW") ~ "Catégorie", # condition ~ résultat
    term == "(Intercept)" ~ "Intercept",  # condition ~ résultat
    TRUE ~ "numérique" # pour toute autre condition (TRUE) ~ résultat
  )) %>% 
  ggplot(mapping = aes(x = estimate, y = term)) +
  geom_vline(xintercept = 0, lty = 2) +
  geom_segment(mapping = aes(x = conf.low, xend = conf.high, yend = term)) +
  geom_point() +
  labs(x = "Coefficient standardisé", y = "") +
  facet_wrap(~type, scales = "free", ncol = 1, strip.position = "right")
```

On y voit qu'à l'exception de la variable `long`, tous les coefficients sont éloignés de 0. Le coefficient `bv` est négatif, indiquant que plus la valeur de `bv` est élevé (donc plus le sol est pauvre en matière organique), plus le rendement est faible. Plus la latitude est élevée (plus on se dirige vers le Nord de l'Argentine), plus le rendement est élevé. La dose d'azote a aussi un effet statistique positif sur le rendement.

Quant aux catégories topographiques, elles sont toutes éloignées de la catégorie `E`, placée à zéro. De plus, les intervalles de confiance à 0.95 ne se chevauchant pas, on peut conclure que la variabilité du phénomène échantillonné n'est pas suffisante pour expliquer les différences importantes d'une à l'autre.

On pourra retrouver des cas où l'effet combiné de plusieurs variables diffère de l'effet des deux variables prises séparément. Par exemple, on pourrait évaluer l'effet de l'azote et celui de la topographie dans un même modèle, puis y ajouter une interaction entre l'azote et la topographie, qui définira des effets supplémentaires de l'azote selon chaque catégorie topographique. C'est ce que l'on appelle une interaction.

Dans l'interface-formule, l’interaction entre l'azote et la topographie est notée `nitro:topo`. Pour ajouter cette interaction, la formule deviendra `yield ~ nitro + topo + nitro:topo`. Une approche équivalente est d'utiliser le raccourci `yield ~ nitro * topo`.

```{r biostats-multivariate-lm-interaction}
modlin_5_sc <- lm(yield ~ nitro * topo,
               data = lasrosas.corn_sc)
summary(modlin_5_sc)
```

Les résultats montre des effets de l'azote et des catégories topographiques, mais il y a davantage d'incertitude sur les interactions, indiquant que l'effet statistique de l'azote est sensiblement le même indépendamment des niveaux topographiques.

Dans le cas des régressions multiples, les résidus ne peuvent pas être présentés selon une variable explicative $x$, puisqu'il y en a plusieurs. On fera l'analyse des résidus selon la variable réponse $y$.

```{r biostats-residuals-mult-reg}
tibble(
  y = lasrosas.corn_sc$yield,
  residus = residuals(modlin_5_sc)
) %>% 
  ggplot(aes(x = y, y = residus)) +
  geom_point() +
  labs(x = "y", y = "Résidus") +
  geom_hline(yintercept = 0, col = "red", size = 1)
```
Dans ce modèle, il y a clairement une structure qui nous échappe! L'ajout de d'autres variables nous permettrait éventuellement d'obtenir une distribution qui s'approche d'un bruit.

#### Les interactions

Une interaction est un effet supplémentaire qui est investigué pour des combinaisons de variables. L'interaction entre l'azote et la topographie est une nouvelle variable créée par la multiplication de l'azote, une variable numérique, et de la topographie, qui ici est une variable catégorielle.

```{r biostats-multivariate-modmat-interaction}
model.matrix(~ nitro * topo, data = lasrosas.corn_sc) %>% head()
```
L'entête de la matrice modèle montre que l'interaction est l'addition de trois variables, qui sont nulles si la catégorie topographique est absente, mais qui prend la dose d'azote pour la catégorie présente seulement.

L'interprétation d'une interaction est spécifique au modèle utilisé. Une manière de l'interpréter est de se demander dans quelles unités elle est exprimée. Dans notre exemple, il s'agit de kg/ha standardisés.

Prenons un autre exemple, cette fois-ci avec des données fictives. Une enquête a été menée, où des personnes évaluait le karma (échelle 0 à 10) de pieds nus, en bas (chaussettes) et/ou en sandales.

```{r biostats-karma-load-data}
karma_df <- read_csv("data/karma_df.csv")
```

Nous désirons savoir quelle est l'effet des bas et des sandales sur le karma, donc 💖 ~ 👡 + 🧦.

```{r biostats-karma-lm-simple}
tidy(lm(karma ~ sandales + bas, karma_df))
```
À partir du scénario à pieds nus d'un karma de 3.67, les sandales ajoutent 2.29 de points de karma, alors que les bas en ajoutent 1.8. Mais ce modèle est incomplet, cas on n'évalue pas l'effet des bas ET des sandales, donc 💖 ~ 👡 * 🧦.

```{r biostats-karma-lm-interaction}
tidy(lm(karma ~ sandales * bas, karma_df))
```

Le modèles est plus clair. Sans interaction, les effets sur le karma des bas et des sandales étaient négativement affectés par l'effet d'interaction `sandales:bas`, le karma étant poussé à la baisse par le [bas blanc dans vos sandales](https://youtu.be/_whvVXX0hCk?t=86).

Il est possible d'ajouter des interactions doubles, triples, quadruples, etc. Mais plus il y a d’interactions, plus votre modèle comprendra de variables et vos tests d'hypothèse perdront en puissance statistique.

#### Les modèles linéaires généralisés

Dans un modèle linéaire ordinaire, un changement constant dans les variables explicatives résulte en un changement constant de la variable-réponse. Cette supposition ne serait pas adéquate si la variable-réponse était un décompte, si elle est booléenne ou si, de manière générale, la variable-réponse ne suivait pas une distribution continue. Ou, de manière plus spécifique, il n'y a pas moyen de retrouver une distribution normale des résidus? On pourra bien sûr transformer les variables (sujet du chapitre 6, en développement). Mais il pourrait s'avérer impossible, ou tout simplement non souhaitable de transformer les variables. Le modèle linéaire généralisé (MLG, ou *generalized linear model* - GLM) est une généralisation du modèle linéaire ordinaire chez qui la variable-réponse peut être caractérisé par une distribution de Poisson, de Bernouilli, etc.

Prenons d'abord cas d'un décompte de vers fil-de-fer (`worms`) retrouvés dans des parcelles sous différents traitements (`trt`). Les décomptes sont typiquement distribué selon une loi de Poisson.

```{r biostats-multivariate-glm-hist}
cochran.wireworms %>% ggplot(aes(x = worms)) + geom_histogram(bins = 10)
```

Explorons les décomptes selon les traitements.

```{r biostats-multivariate-glm-boxplot}
cochran.wireworms %>% ggplot(aes(x = trt, y = worms)) + geom_boxplot()
```

Les traitements semble à première vue avoir un effet comparativement au contrôle. Lançons un MLG avec la fonction `glm()`, et spécifions que la sortie est une distribution de Poisson. Bien que la fonction de lien (`link = "log"`) soit explictement imposée, le log est la valeur par défaut pour les distributions de Poisson. Ainsi, les coefficients du modèles devront être interprétés selon un modèle $log \left(worms \right) = intercept + pente \times coefficient$.

```{r biostats-multivariate-glm-model}
modglm_1 <- glm(worms ~ trt, cochran.wireworms, family = stats::poisson(link="log"))
summary(modglm_1)
```

L'interprétation spécifique des coefficients d'une régression de Poisson doit passer par la fonction de lien $log \left(worms \right) = intercept + pente \times coefficient$. Le traitement de référence (K), qui correspond à l'intercept, sera accompagné d'un nombre de vers de $exp \left(0.1823\right) = 1.20$ vers, et le traitement M, à $exp \left(1.6422\right) = 5.17$ vers. Cela correspond à ce que l'on observe sur les boxplots plus haut.

Il est très probable (p-value de ~0.66) qu'un intercept (traitement K) de 0.18 ayant une erreur standard de 0.4082 ait été généré depuis une population dont l'intercept est nul. Quant aux autres traitements, leurs effets sont tous significatifs au seuil 0.05, mais peuvent-ils être considérés comme équivalents?

```{r biostats-multivariate-glm-ci}
intervals <- tibble(Estimate = coefficients(modglm_1), # [-1] enlever l'intercept
                    LL = confint(modglm_1)[, 1], # [-1, ] enlever la première ligne, celle de l'intercept
                    UL = confint(modglm_1)[, 2],
                    variable = names(coefficients(modglm_1)))
intervals
```

```{r biostats-multivariate-glm-plot-effects}
ggplot(data = intervals, mapping = aes(x = Estimate, y = variable)) +
    geom_vline(xintercept = 0, lty = 2) +
    geom_segment(mapping = aes(x = LL, xend = UL,
                               y = variable, yend = variable)) +
    geom_point() +
    labs(x = "Coefficient", y = "")
```

Les intervalles de confiance se superposant, on ne peut pas conclure qu'un traitement est lié à une réduction plus importante de vers qu'un autre, au seuil 0.05.

Maintenant, à défaut de trouver un tableau de données plus approprié, prenons le tableau `mtcars`, qui rassemble des données sur des modèles de voitures. La colonne `vs`, pour v-shaped, inscrit 0 si les pistons sont droit et 1 s'ils sont placés en V dans le moteur. Peut-on expliquer la forme des pistons selon le poids du véhicule (`wt`)?

```{r biostats-multivariate-glm-cars}
mtcars %>% sample_n(6)
```

```{r biostats-multivariate-glm-cars-plot}
mtcars %>%
    ggplot(aes(x = wt, y = vs)) + geom_point()
```

Il semble y avoir une tendance: les véhicules plus lourds ont plutôt des pistons droits (`vs = 0`). Vérifions cela.

```{r biostats-multivariate-glm-cars-model}
modglm_2 <- glm(vs ~ wt, data = mtcars, family = stats::binomial())
summary(modglm_2)
```

**Exercice**. Analyser les résultats.

#### Les modèles non-linéaires

La hauteur d'un arbre en fonction du temps n'est typiquement pas linéaire. Elle tend à croître de plus en plus lentement jusqu'à un plateau. De même, le rendement d'une culture traité avec des doses croissantes de fertilisants tend à atteindre un maximum, puis à se stabiliser.

Ces phénomènes ne peuvent pas être approximés par des modèles linéaires. Examinons les données du tableau `engelstad.nitro`.

```{r biostats-multivariate-nls-data}
engelstad.nitro %>% sample_n(10)
```

```{r biostats-multivariate-nls-plot}
engelstad.nitro %>%
    ggplot(aes(x = nitro, y = yield)) +
        facet_grid(year ~ loc) +
        geom_line() +
        geom_point()
```

Le modèle de Mitscherlich pourrait être utilisé.

$$ y = A \left( 1 - e^{-R \left( E + x \right)} \right) $$

où $y$ est le rendement, $x$ est la dose, $A$ est l'asymptote vers laquelle la courbe converge à dose croissante, $E$ est l'équivalent de dose fourni par l'environnement et $R$ est le taux de réponse.

Explorons la fonction.

```{r biostats-multivariate-nls-miscth-explore}
mitscherlich_f <- function(x, A, E, R) {
    A * (1 - exp(-R*(E + x)))
}

x <- seq(0, 350, by = 5)
y <- mitscherlich_f(x, A = 75, E = 30, R = 0.02)

ggplot(tibble(x, y), aes(x, y)) +
    geom_point(data = engelstad.nitro, aes(x = nitro, y = yield)) +
    geom_line() + ylim(c(0, 100))
```

**Exercice**. Changez les paramètres pour visualiser comment la courbe réagit.

Nous pouvons décrire le modèle grâce à l'interface formule dans la fonction `nls()`. Notez que les modèles non-linéaires demandent des stratégies de calcul différentes de celles des modèles linéaires. En tout temps, nous devons identifier des valeurs de départ raisonnables pour les paramètres dans l'argument `start`. Vous réussirez rarement à obtenir une convergence du premier coup avec vos paramètres de départ. Le défi est d'en trouver qui permettront au modèle de converger. Parfois, le modèle ne convergera jamais. D'autres fois, il convergera vers des solutions différentes selon les variables de départ choisies.

```
modnl_1 <- nls(yield ~ A * (1 - exp(-R*(E + nitro))),
                data = engelstad.nitro,
                start = list(A = 50, E = 10, R = 0.2))
```

Le modèle ne converge pas (le bloc de calcul est désactivé). Essayons les valeurs prises plus haut, lors de la création du graphique, qui semblent bien s'ajuster.

```{r biostats-multivariate-nls-model}
modnl_1 <-  nls(yield ~ A * (1 - exp(-R*(E + nitro))),
                data = engelstad.nitro,
                start = list(A = 75, E = 30, R = 0.02))
```

Bingo! Voyons maintenant le sommaire.

```{r biostats-multivariate-nls-summary}
summary(modnl_1)
```

Les paramètres sont différents de zéro, et donnent la courbe suivante.

```{r biostats-multivariate-nls-resultat-plot}
x <- seq(0, 350, by = 5)
y <- mitscherlich_f(x,
                    A = coefficients(modnl_1)[1],
                    E = coefficients(modnl_1)[2],
                    R = coefficients(modnl_1)[3])

ggplot(tibble(x, y), aes(x, y)) +
    geom_point(data = engelstad.nitro, aes(x = nitro, y = yield)) +
    geom_line() + ylim(c(0, 100))
```

Et les résidus...

```{r biostats-multivariate-nls-residuals}
tibble(res = residuals(modnl_1)) %>%
    ggplot(aes(x = res)) + geom_histogram(bins = 20)
```

```{r biostats-multivariate-nls-residuals2}
tibble(nitro = engelstad.nitro$nitro, res = residuals(modnl_1)) %>%
    ggplot(aes(x = nitro, y = res)) +
        geom_point() +
        geom_hline(yintercept = 0, colour = "red")
```

Les résidus ne sont pas distribués normalement, mais semble bien partagés de part et d'autre de la courbe.

### Modèles à effets mixtes

Lorsque l'on combine des variables fixes (testées lors de l'expérience) et des variables aléatoire (variation des unités expérimentales), on obtient un modèle mixte. Les modèles mixtes peuvent être univariés, multivariés, linéaires ordinaires ou généralisés ou non linéaires.

À la différence d'un effet fixe, un effet aléatoire sera toujours distribué normalement avec une moyenne de 0 et une certaine variance. Dans un modèle linéaire où l'effet aléatoire est un décalage d'intercept, cet effet s'additionne aux effets fixes:

$$ y = X \beta + Z b + \epsilon $$

où:

$Z$ est la matrice du modèle à $n$ observations et $p$ variables aléatoires. Les variables aléatoires sont souvent des variables nominales qui subissent un encodage catégoriel.

$$ Z = \left( \begin{matrix}
z_{11} & \cdots & z_{1p}  \\
z_{21} & \cdots & z_{2p}  \\
\vdots & \ddots & \vdots  \\
z_{n1} & \cdots & z_{np}
\end{matrix} \right) $$

$b$ est la matrice des $p$ coefficients aléatoires.

$$ b = \left( \begin{matrix}
b_0  \\
b_1  \\
\vdots \\
b_p
\end{matrix} \right) $$

Le tableau `lasrosas.corn`, utilisé précédemment, contenait trois répétitions effectués au cours de deux années, 1999 et 2001. Étant donné que la répétition R1 de 1999 n'a rien à voir avec la répétition R1 de 2001, on dit qu'elle est **emboîtée** dans l'année.

Le module `nlme` nous aidera à monter notre modèle mixte.

```{r biostats-multivariate-nlme-model}
library("nlme")
mmodlin_1 <- lme(fixed = yield ~ lat + long + nitro + topo + bv,
                 random = ~ 1|year/rep,
                 data = lasrosas.corn)
```

À ce stade vous devriez commencer à être familier avec l'interface formule et vous deviez saisir l'argument `fixed`, qui désigne l'effet fixe. L'effet aléatoire, `random`, suit un tilde `~`. À gauche de la barre verticale `|`, on place les variables désignant les effets aléatoire sur la pente. Nous n'avons pas couvert cet aspect, alors nous le laissons à `1`. À droite, on retrouve un structure d'emboîtement désignant l'effet aléatoire: le premier niveau est l'année, dans laquelle est emboîtée la répétition.

```{r biostats-multivariate-nlme-model}-summ
summary(mmodlin_1)
```

La sortie est semblable à celle de la fonction `lm()`.

#### Modèles mixtes non-linéaires

Le modèle non linéaire créé plus haut liait le rendement à la dose d'azote. Toutefois, les unités expérimentales (le site `loc` et l'année `year`) n'étaient pas pris en considération. Nous allons maintenant les considérer.

Nous devons décider la structure de l'effet aléatoire, et sur quelles variables il doit être appliqué - la décision appartient à l'analyste. Il me semble plus convenable de supposer que le site et l'année affectera le rendement maximum plutôt que l'environnement et le taux: les effets aléatoires seront donc affectés à la variable `A`. Les effets aléatoires n'ont pas de structure d'emboîtement. L'effet de l'année sur A sera celui d'une pente et l'effet de site sera celui de l'intercept. La fonction que nous utiliserons est `nlme()`.

```{r biostats-multivariate-nlme-nl-model}
mm <- nlme(yield ~ A * (1 - exp(-R*(E + nitro))),
           data = engelstad.nitro,
           start = c(A = 75, E = 30, R = 0.02),
           fixed = list(A ~ 1, E ~ 1, R ~ 1),
           random = A ~ year | loc)
summary(mm)
```

Et sur graphique:

```{r biostats-multivariate-nlme-nl-plot}
engelstad.nitro %>%
  ggplot(aes(x = nitro, y = yield)) +
  facet_grid(year ~ loc) +
  geom_line(data = tibble(nitro = engelstad.nitro$nitro,
                          yield = predict(mm, level = 0)),
            colour = "grey35") +
  geom_point() +
  ylim(c(0, 95))
```


Les modèles mixtes non linéaires peuvent devenir très complexes lorsque les paramètres, par exemple A, E et R, sont eux-même affectés linéairement par des variables (par exemple `A ~ topo`). Pour aller plus loin, consultez [Parent et al. (2017) ](https://doi.org/10.3389/fenvs.2017.00081) ainsi que les [calculs associés à l'article](https://github.com/essicolo/site-specific-multilevel-modeling-of-potato-response-to-nitrogen-fertilization). Ou écrivez-moi un courriel pour en discuter!

**Note**. L'interprétation de p-values sur les modèles mixtes est controversée. À ce sujet, Douglas Bates a écrit une longue lettre à la communauté de développement du module `lme4`, une alternative à `nlme`, qui remet en cause l'utilisation des p-values, [ici](https://stat.ethz.ch/pipermail/r-help/2006-May/094765.html). De plus en plus, pour les modèles mixtes, on se tourne vers les statistiques bayésiennes, couvertes dans le chapitre \@ref(chapitre-biostats-bayes) avec le module greta. Mais en ce qui a trait aux modèles mixtes, le module [`brms`](https://github.com/paul-buerkner/brms) automatise bien des aspects de l'approche bayésienne.

### Aller plus loin

#### Statistiques générales:
- [The analysis of biological data](https://ariane25.bibl.ulaval.ca/ariane/wicket/detail?c=ariane&m=S&rq.ct=PE&rq.fa=false&rq.r.esc=false&rq.r.l%5B0%5D.c=TI&rq.r.l%5B0%5D.ex=false&rq.r.l%5B0%5D.op=AND&rq.r.l%5B0%5D.v=the+analysis+of+biological+data&rq.r.la=*&rq.r.loc=*&rq.r.pft=true&rq.r.ta=*&rq.r.td=*&rq.rows=5&rq.st=1)

#### Statistiques avec R

- Disponibles en version électronique à la bibliothèque de l'Université Laval:
    - Introduction aux statistiques avec R: [Introductory statistics with R](https://ariane25.bibl.ulaval.ca/ariane/wicket/detail?c=ariane&m=S&rq.ct=PE&rq.fa=false&rq.r.esc=false&rq.r.l%5B0%5D.c=TI&rq.r.l%5B0%5D.ex=false&rq.r.l%5B0%5D.op=AND&rq.r.l%5B0%5D.v=Introductory+statistics+with+R&rq.r.la=*&rq.r.loc=*&rq.r.pft=true&rq.r.ta=*&rq.r.td=*&rq.rows=1&rq.st=0)
    - Approfondir les statistiques avec R: [The R Book, Second edition](https://ariane25.bibl.ulaval.ca/ariane/wicket/detail?c=ariane&m=S&rq.ct=PE&rq.fa=false&rq.r.esc=false&rq.r.l%5B0%5D.c=TI&rq.r.l%5B0%5D.ex=false&rq.r.l%5B0%5D.op=AND&rq.r.l%5B0%5D.v=the+r+book&rq.r.la=*&rq.r.loc=*&rq.r.pft=true&rq.r.ta=*&rq.r.td=*&rq.rows=15&rq.st=2)
    - Approfondir les modèles à effets mixtes avec R: [Mixed Effects Models and Extensions in Ecology with R](https://ariane25.bibl.ulaval.ca/ariane/wicket/detail?c=ariane&m=S&rq.ct=PE&rq.fa=false&rq.r.esc=false&rq.r.l%5B0%5D.c=TI&rq.r.l%5B0%5D.ex=false&rq.r.l%5B0%5D.op=AND&rq.r.l%5B0%5D.v=Mixed+Effects+Models+and+Extensions+in+Ecology+with+R&rq.r.la=*&rq.r.loc=*&rq.r.pft=false&rq.r.ta=*&rq.r.td=*&rq.rows=2&rq.st=1)
- [ModernDive](https://moderndive.com/index.html), un livre en ligne offrant une approche moderne avec le package `moderndive`.

```{r biostats-rm-all, include=FALSE}
rm(list = ls())
```