Atelier 7: Modèles linéaires et généralisés linéaires mixtes

.title[
# Atelier 7: Modèles linéaires et généralisés linéaires mixtes
]
.subtitle[
## Série d’ateliers R du CSBQ
]
.author[
### Centre de la Science de la Biodiversité du Québec
]

---

# À propos de cet atelier

[![badge](https://img.shields.io/static/v1?style=flat-square&label=Diapos&message=07&color=red&logo=html5)](https://r.qcbs.ca/workshop07/pres-fr/workshop07-pres-fr.html) 
[![badge](https://img.shields.io/static/v1?style=flat-square&label=livre&message=07&logo=github)](https://r.qcbs.ca/workshop07/book-fr/index.html) 
[![badge](https://img.shields.io/static/v1?style=flat-square&label=script&message=07&color=2a50b8&logo=r)](https://r.qcbs.ca/workshop07/book-fr/workshop07-script-fr.R) 
[![badge](https://img.shields.io/static/v1?style=flat-square&label=repo&message=dev&color=6f42c1&logo=github)](https://github.com/QCBSRworkshops/workshop07)

???
Note: Rappeler aux participants que le matériel de l'atelier est disponible sur le site web des ateliers R et qu'il y a un 'Bookdown' avec du texte supplémentaire. Le contenu de l'atelier est plutôt lourd, donc ils pourront retourner aux notes et les diapositives de l'atelier pour revisiter les concepts qu'ils auront moins bien compris.

---

en modifiant et améliorant le matériel dans le cadre des *P*rix d’*A*pprentissage et de *D*éveloppement
]

**2022** - **2021** - **2020**

[Maxime Fraser Franco]()
 
[Hassen Allegue]()
 
[Linley Sherin]()
 
[Pedro Henrique P. Braga]()
 
[Katherine Hébert]()
 
[Kevin Cazelles]()
 
[Janaína Serrano]()
 
[Dominique Caron]()
]
]

[Nicolas Princeloup]()
 
[Marie Hélène Brice]()

**2016** - **2015** - **2014**

[Catherine Baltazar]()
 
[Dalal Hanna]()
 
[Jacob Ziegler]()
 
[Cédric Frenette Dussault]()
 
[Vincent Fugère]()
 
[Thomas Lamy]()
 
[Zofia Taranu]()

]
]

---

# Matériel requis

Pour suivre cet atelier, il est nécessaire d'avoir téléchargé et installé les dernières versions de [RStudio](https://rstudio.com/products/rstudio/download/#download) et de [R](https://cran.rstudio.com/).

* [lme4](http://cran.r-project.org/web/packages/lme4/index.html)
* [MASS](http://cran.r-project.org/web/packages/MASS/index.html)
* [vcdExtra](http://cran.r-project.org/web/packages/vcdExtra/index.html)
* [bbmle](http://cran.r-project.org/web/packages/bbmle/index.html)
* [MuMIn](http://cran.r-project.org/web/packages/MuMIn/index.html)
* [ggplot2](http://cran.r-project.org/web/packages/ggplot2/index.html)
* [DescTools](http://cran.r-project.org/web/packages/DescTools/index.html)
* [remotes](http://cran.r-project.org/web/packages/remotes/index.html)
* [gridExtra](http://cran.r-project.org/web/packages/gridExtra/index.html)
* [lattice](http://cran.r-project.org/web/packages/lattice/index.html)

]

Pour les installer à partir du CRAN, exécutez :

```r
install.packages(c("lme4",
                      "MASS",
                      "vcdExtra",
                      "bbmle", 
                      "MuMIn",
                      "ggplot2",
                      "DescTools",
                      "remotes",
                      "gridExtra",
                      "lattice"))
```
]

???
Note: Pour sauver du temps: (1) Demander aux participants d'installer les paquets avant l'atelier. (2) Ne pas attendre que tout le monde réussisse à les installer, et aider ceux qui ont rencontrer des problèmes lors du défi #1.

---

# Matériel requis

- <a href="https://r.qcbs.ca/workshop07/pres-en/data/fishdata.csv">fishdata.csv</a> 
- <a href="https://r.qcbs.ca/workshop07/pres-en/data/arabidopsis.csv">arabidopsis.csv</a>
- <a href="https://r.qcbs.ca/workshop07/pres-en/data/inverts.csv">inverts.csv</a>
]

.pull-left2[
Tout au long de cet atelier, il y aura une série de **défis** que vous pouvez reconnaître par ce cube rubix.
]

]

---

# Objectifs d'apprentissage

**1.** Décrire les modèles (généralisés) à effets mixtes

**2.** Identifier les situations dans lesquelles l'utilisation d'effets mixtes est appropriée

**3.** Mettre en œuvre des modèles linéaires mixtes de base avec `R`

**4.** Exécuter des modèles linéaires généralisés mixtes de base avec `R`

**5.** Valider, interpréter et visualiser les modèles mixtes avec `R`

---
# Question de recherche

Pour répondre, nous utiliserons un jeu de données où la longueur corporelle de 3 espèces de poissons (10 individus par espèce) a été mesurée dans 6 lacs différents.

---
# Défi 1 ![:cube]()

**Familiarisez-vous avec le jeu de données**

**1.** Ouvrez le script de l'atelier dans `R`

**2.** Ouvrez le jeu de données dans `R`

**3.** Reproduisez les graphiques 1 à 3 (dans le script) de la relation entre la position trophique et la taille. Observez les graphiques, puis évaluez ce que vous observez.

???
Note: Ce défi devrait être assez rapide. Peut-être simplement accorder un peu de temps de réflexion individuel pendant que vous aidez les participants qui ont eu de la difficulté à installer les paquets. Si il y en a aucun, vous pouvez simplement avoir une discussion de groupe un graphique à la fois.

---
exclude: true

# Familiarisez-vous avec le jeu de données

**1.** Ouvrez le jeu de données dans `R`

**2.** Ouvrez le script de l'atelier dans `R`

**3.** Visualisez la relation entre taille et position trophique

---
# Solution ![:cube]()

---
# Solution ![:cube]()

---
# Solution ![:cube]()

---
# Discussion de groupe ![:cube]()

**Est-ce qu'on s'attend à ce que la position trophique augmente avec la longueur corporelle exactement de la même façon pour : **

* toutes les espèces?
* tous les lacs?

**Comment ces relations pourraient-elles différer?**

---
# Pourquoi choisir un MLM?

### Les données biologiques et écologiques peuvent être complexes!

* Peuvent contenir une structure hiérarchique

* Plusieurs covariables et facteurs

* Échantillons/design expérimental non équilibrés

???
Note: Utiliser la base de données comme example. Nous avons des données de longueur de poissons et de niveau toophique pour 3 espèces et 6 lacs: cela créer une structure hiérarchique où nous avons un niveau divisant les données par espèces et un autre par lac.

---
# Pourquoi choisir un MLM?

### Comment pourrions-nous analyser ces données?

**Option 1. Séparer**

- Faire une analyse séparée pour chaque espèce dans chaque lac

**Option 2. Regrouper**

- Faire une seule analyse en ignorant les variables espèce et lac

---
# Pourquoi choisir un MLM?

On pourrait faire une analyse pour chaque espèce :
* Estimer 6 ordonnées à l'origine et 6 pentes pour chaque espèce (i.e. 6 lacs)

* Taille d'échantillon *n* = 10 pour chaque analyse (i.e. 10 poissons/espèce/lac)

* Peu de chances de détecter un effet en raison de la faible taille d’échantillon *n*
]

---
# Pourquoi choisir un MLM?

* Et la pseudoreplication? (les poissons d'un même lac et d'une même espèce pourraient être corrélés)

* Il y a beaucoup de bruit dans les données! Une partie pourrait être expliquée par les **différences** entre **les espèces** et les **lacs**.
]

---
# Pourquoi choisir un MLM?

Pour **notre question**, on veut seulement savoir s'il y a un **effet général de la longueur corporelle sur la position trophique**.

Par contre :

* **Cette relation pourrait varier** par **espèce** à cause de différents processus biologiques n'ayant pas été mesurés (ex: taux de croissance) et/ou par **lac** à cause de variables environnementales non-mesurées (ex: différences dans la disponibilité de nourriture).

* **On doit contrôler** ces effets potentiels dans le modèle.

Les MLMs permettent de **séparer et regrouper** l'analyse. Ils:

1. Prennent en compte la variabilité spécifique à chaque espèce et chaque lac (**séparer**) tout en calculant moins de paramètres qu'une régression classique.

2. Utilisent toutes les données disponibles (**regrouper**) tout en contrôlant les différences entre les lacs et les espèces (pseudo-réplication).

---
# Pourquoi choisir un MLM?

### Effet fixe vs effet aléatoire

Dans la littérature des MLMs, vous rencontrerez souvent ces termes.

Il existe plusieurs façons de les décrire et nous vous présenterons ici celles que nous trouvons les plus faciles à appliquer.

---
# Pourquoi choisir un MLM?

### Effet fixe : processus déterministes

Les données proviennent : 
* de tous les niveaux possibles d'un facteur (**variable qualitative**) 
* d'un prédicteur (**variable quantitative**)

On souhaite émettre des conclusions à propos des niveaux du facteur ou de la relation entre le prédicteur et la variable réponse.

*Exemple : Nous sommes intéressés par l'effet d'un traitement (ex: coupes à blanc) sur une variable réponse (ex: richesse spécifique). La variable traitement est le facteur à deux niveaux: forêt coupée à blanc ou non coupée.*

???
Note: Mettre l'emphase sur le fait que nous avons échantillonné tous les niveaux d'intérêts du facteurs. Nous ne sommes pas intéressés par les autres niveaux non échantillonnés.

---
# Pourquoi choisir un MLM?

### Effet aléatoire : processus stochastiques

* Seulement des **variables qualitatives** = facteur aléatoire

* Les données ne comprennent qu'un échantillon aléatoire de tous les niveaux de facteurs possibles, tous d'intérêt

* Permet de structurer le processus d'erreur

*Exemple : Nous échantillons des poissons dans 6 de 15 lacs d'une région. La variable lac est le facteur avec 15 niveaux. Bien que nous n'avons qu'échantillonné 6 niveaux, nous voulons pouvoir faire des conclusions sur la variabilité entre tous les lacs de la région.*

???
Note: Mettre l'emphase sur le fait que nous n'avons pas nécessairement des données pour tous les niveaux d'intérêt du facteur aléatoire.

---
# Pourquoi choisir un MLM?

### Comment fonctionnent les MLMs?

**A.** Les ordonnées à l'origine et/ou les pentes peuvent être considérées comme propres (c.-à.-d. variant en fonction de) à une population donnée (**effet aléatoire**) (ex: par lac et/ou par espèce).

**B.** Les ordonnées à l'origine, les pentes et leur intervalle de confiance sont ajustés pour **prendre en compte la structure des données**.

---
# Effet aléatoire : faire varier les ordonnées à l'origine et/ou les pentes

* On suppose que les ordonnées à l'origine et les pentes proviennent d'une distribution normale.
* Le modèle estime une ordonnée à l'origine et une pente moyenne ainsi que l'écart type de l'effet aléatoire (distribution normale).

Évite d'estimer une ordonnée à l'origine et une pente par espèce (+2 paramètres par espèce) = **économise des degrés de liberté** (moins de paramètres à estimer).

???
Note: Cela est assez important. Passez du temps pour bien expliquer the l'on assume the l'ordonnée et/ou la pente de chaque niveau proviennent d'une même distribution (ils ne sont pas totalement indépendent), ce qui permet de réduire le nombre de paramètres à estimer par le modèle.

---
# Effet aléatoire sur l'ordonnée à l'origine

On fait la supposition que les ordonnées à l'origine proviennent d'une distribution normale.

On a seulement besoin d'estimer la moyenne (ordonnée à l'origine générale) et l'écart type de la distribution normale (effet aléatoire) au lieu d'estimer une ordonnée à l'origine par espèce (ce qui ajouterait 2 paramètres).

.comment[Notez que plus votre facteur a de niveaux, plus la moyenne et l'écart-type de la distribution normale seront estimés avec précision. Trois niveaux est un peu faible, mais c'est certainement plus facile à visualiser !]

???
J'ai enlevé ceci pour rajouter la note présente dans la présentation en anglais.
Au lieu d'estimer une ordonnée à l'origine par espèce (3 paramètres), on en estime une générale ainsi qu'un effet aléatoire (2 paramètres). Avec `$n$` espèce, dans le premier cas on estime `$n-1$` paramètres alors qu'on reste à 2 paramètres dans le second cas!

---
# Effet aléatoire sur l'ordonnée à l'origine

<img style="float: right; width:50%;margin: 1%" src="images/fig_8_w5.png">
Même principe pour les lacs.

Estime 2 paramètres (moyenne et écart-type) au lieu de 6 ordonnées à l'origine.

Cela **économise des degrés de liberté** (moins de paramètres à estimer).

---
# Effet aléatoire sur la pente

<img style="float: right; width:50%;margin: 1%" src="images/fig_9_w5.png">
Le même principe s’applique aux pentes qui varient selon un facteur donné.

Comme pour les ordonnées à l'origine, seuls la moyenne et l’écart-type des pentes sont estimés au lieu de trois pentes distinctes.

---
# Tenir compte de la structure des données
 
**Qu'arrive-t-il si j'ai peu d'échantillons (faible `$n$`) pour un niveau des facteurs?**

Si une espèce ou un lac est peu représenté dans les données, le modèle va accorder plus d'importance au modèle groupé pour estimer l'ordonnée à l'origine et la pente de cette espèce ou de ce lac (processus de « shrinkage »).

Idéalement, il faudrait avoir un minimum de `$n$`=3 par niveau d'un facteur.

---
# Tenir compte de la structure des données
 
**Comment évaluer l'impact d'un effet aléatoire sur le modèle?**

* Les intervalles de confiance des ordonnées à l'origine et des pentes **générales** sont ajustés pour tenir compte de la pseudo-réplication basée sur le **coefficient de corrélation intra-classe (CCI)**.

* Le CCI permet de savoir la proportion de variation dans la variable réponse qui est expliquée par l'effet aléatoire, tel que :

`$$CCI = \frac{\sigma_{\alpha}^2}{\sigma_{\alpha}^2 + \sigma_{\varepsilon}^2}$$`
Le CCI correspond au **ratio** entre la variance d'un effet aléatoire (e.g. ordonnées à l'origine des espèces) et la variance totale.

Le CCI nous indique donc **à quel point** la position trophique moyenne entre chaque espèce ou chaque lac (c.-à-d. les ordonnées à l'origine) **varie**.

???

Note: Dire que différentes notations sont possibles selon le livre/article et comment l'équation du modèle est écrite.

---
# Tenir compte de la structure des données

![:scale 90%](images/CCIplot2_qcbs.png)

Le % de variance (CCI) est élevé, car **les espèces diffèrent fortement** dans leur position trophique moyenne.

Les intervalles de confiance pour la pente et l'ordonnée à l'origine générale sont grands.
]

![:scale 85%](images/CCIplot_qcbs.png)

Le % de variance (CCI) est faible, car **les espèces diffèrent peu** dans leur position trophique moyenne.

Les intervalles de confiance pour la pente et l'ordonnée à l'origine générale sont petits.
]

---
exclude: true
# Tenir compte de la structure des données

![](images/fig_10_w5.png)

Les points provenant d'un même lac sont traités comme une seule observation car **très corrélés**.

![:faic](arrow-right) petite taille effective de l'échantillon et grands intervalles de confiance pour la pente et l'ordonnée à l'origine.
]

![](images/fig_11_w5.png)

Les points provenant d'un même lac sont traités indépendamment car **peu corrélés**.

![:faic](arrow-right) grande taille effective de l'échantillon et petits intervalles de confiance pour la pente et l'ordonnée à l'origine.
]

---
# Défi 2 ![:cube]()

**Comment le CCI et l'intervalle de confiance seront-ils affectés dans ces deux scénarios?**

**Q1.** Les positions trophiques des poissons ne varient pas entre les lacs.

**Q2.** Les positions trophiques des poissons sont similaires au sein d'un même lac, mais différentes entre les lacs.

???

Note: Pour ce défi, vous pouvez créer un sondage Zoom. Cela ne devrait pas être trop long.

---
# Solution ![:cube]()

**Q1.** Les positions trophiques des poissons ne varient pas entre les lacs.

**Q2.** Les positions trophiques des poissons sont similaires au sein d'un même lac, mais différentes entre les lacs.

**Pour plus de détails sur le CCI** : 
Nakagawa et Schielzeth (2013)
https://doi.org/10.1111/j.1469-185X.2010.00141.x
Nakagawa *et al.* (2017) 
https://doi.org/10.1098/rsif.2017.0213

???

Note: Répéter que le CCI mesure le ratio de la variance "expliquée" par le facteur aléatoire. De plus, plus les intervalles de confiance pour l'ordonnée à l'origine et la pente globale est grande, plus la variation interclasse est grande.

---
# Comment implémenter un MLM avec R?

###### **Étape 1.** Construction du modèle *a priori* et exploration des données

**Étape 2.** Coder les modèles potentiels et sélection du meilleur modèle

**Étape 3.** Validation du modèle

**Étape 4.** Interprétation et visualisation des résultats

---
# Étape 1. construction du modèle *a priori*

**Modèle basé sur connaissance *a priori*:**

* Nous voulons déterminer si la position trophique peut être prédite par la longueur corporelle, tout en prenant en compte la variation entre les espèces et les lacs.

* Nous voulons donc un modèle qui ressemble à ceci :

`$$PT_{ijk} \sim Longueur_i + Lac_j + Espèce_k + \varepsilon_{ijk}$$`

???

Note: Décrire chaque élément de l'équation. La position trophique est prédite par la longueur, le lac de provenance et l'espèce du poisson. Le dernier terme décrit l'erreur entre la position trophique prédites par les variables explicatives et la vraie position trophique.

---
# Étape 1. exploration des données

**Les données ont-elles la bonne structure?**

```r
fish.data <- read.csv('data/fishdata.csv')
str(fish.data)
# 'data.frame':	180 obs. of 4 variables:
# $ Lake : chr "L1" "L1" "L1" "L1" ...
# $ Fish_Species: chr "S1" "S1" "S1" "S1" ...
# $ Fish_Length : num 105 195 294 414 237 ...
# $ Trophic_Pos : num 2.6 2.7 2.74 2.74 2.79 ...
```

Il est recommandé de faire le ménage de votre espace de travail (`rm.list()`) avant de construire un modèle.

???

Note: Vous pouvez aller relativement rapidement pour l'étape 1 (diapositives 37-50). Les participants devraient être relativement à l'aise avec ces étapes.

---
# Étape 1. exploration des données

**Regardez la distribution des échantillons pour chaque facteur:**

```r
table(fish.data[ , c("Lake", "Fish_Species")])
#     Fish_Species
# Lake S1 S2 S3
#   L1 10 10 10
#   L2 10 10 10
#   L3 10 10 10
#   L4 10 10 10
#   L5 10 10 10
#   L6 10 10 10
```

Ce jeu de données est parfaitement équilibré, mais les **modèles mixtes peuvent analyser des designs expérimentaux non équilibrés**, comme c'est souvent le cas en écologie!

???

Note: Expliquer ce qu'est un jeu de données équilibré vs non équilibrés (c-a-d. le nombre d'échantillons dans chaque classe.)

---
# Étape 1. exploration des données

**Regardez la distribution des variables continues :**

```r
hist(fish.data$Fish_Length, xlab = "Longueur (mm)", main = "")
hist(fish.data$Trophic_Pos, xlab = "Position trophique", main = "")
```

Des déviations majeures pourraient causer des problèmes d'hétéroscédasticité. Si nécessaire, faites des transformations. Dans ce cas-ci, **les données semblent correctes**.

---
# Étape 1. exploration des données

**Vérifiez la colinéarité entre vos variables explicatives**

Le problème avec les prédicteurs colinéaires est simplement qu'ils expliquent la même chose, alors leur effet sur la variable réponse sera confondu dans le modèle.

Dans cet exemple, il n’y a pas de risque de colinéarité, car il y a seulement une variable continue. Si vous aviez une autre variable continue (`Var2`), une façon simple de vérifier la colinéarité serait :

```r
plot(data)

cor(var1, var2)
```

Voici un [exemple de colinéarité](https://yetanotheriteration.netlify.app/2018/01/high-collinearity-effect-in-regressions/).

---
# Défi 3 ![:cube]()

Quelles mesures supplémentaires aurions-nous pu prendre sur le terrain et qui auraient pu être fortement corrélées avec la longueur corporelle?

> Un exemple est la masse du poisson – c’est une variable fortement corrélée avec la longueur du poisson. Par conséquent, nous ne voulons pas inclure ces deux variables dans le même modèle.

???

Note: Ce défi est principalement pour interagir un peu avec les participants. Poser la question à l'audience, et ça ne devrait pas être trop long.

---
# Étape 1. exploration des données

**Considérez l'échelle de vos données**

* Si deux variables dans un même modèle ont des échelles très différentes, il est probable que le modèle indique un problème de convergence en essayant de calculer les paramètres.

* La <a href="https://fr.wikipedia.org/wiki/Cote_Z_(statistiques)">correction Z</a> standardise les variables et résout ce problème (fonction `scale()` dans `R`) :

`$$z = \frac{x-moyenne(x)}{écart.type(x)}$$`
---
# Étape 1. exploration des données

**Considérez l'échelle de vos données**

* Longueur corporelle ![:faic](arrow-right) Longue échelle

* Position trophique ![:faic](arrow-right) Courte échelle

---
# Étape 1. exploration des données

**Considérez l'échelle de vos données**

Étant donné que nos données ont des échelles très différentes, on applique la **correction Z**.

```r
# Longueur corrigée, "à la main"
fish.data$Z_Length <- (fish.data$Fish_Length - mean(fish.data$Fish_Length)) / 
 sd(fish.data$Fish_Length)

# Position trophique corrigée, avec la fonction scale
fish.data$Z_TP <- scale(fish.data$Trophic_Pos)
```

---
# Étape 1. exploration des données

Pour savoir si un modèle mixte est nécessaire pour vos données, vous devez déterminer s'il est important de prendre en compte l'effet aléatoire de facteurs qui pourraient influencer la relation qui vous intéresse (dans notre cas, lac et espèce).

**Nous pouvons le faire en :**

1. Créant un modèle linéaire sans les facteurs qui pourraient avoir un effet aléatoire

2. Calculant les résidus de ce modèle linéaire

3. Produisant un graphique de la valeur des résidus en fonction des niveaux des facteurs potentiellement aléatoires

---
# Étape 1. exploration des données

Créez un modèle linéaire sans les facteurs

```r
lm.test <- lm(Z_TP ~ Z_Length, data = fish.data)
```

Calculez les résidus de ce modèle linéaire

```r
lm.test.resid <- rstandard(lm.test)
```

---
# Étape 1. exploration des données

Représentez graphiquement la valeur des résidus en fonction des niveaux des facteurs

```r
plot(lm.test.resid ~ as.factor(fish.data$Fish_Species),
     xlab = "Species", ylab = "Standardized residuals")

abline(0, 0, lty = 2)

plot(lm.test.resid ~ as.factor(fish.data$Lake),
     xlab = "Lake", ylab = "Standardized residuals")

abline(0, 0, lty = 2)
```

---
# Étape 1. exploration des données

Représentez graphiquement la valeur des résidus en fonction des niveaux des facteurs

.alert[Ces patrons suggèrent qu'il y a de la variance résiduelle qui pourrait être expliquée par ces facteurs. Ils devraient donc être inclus dans le modèle.]

---
# Comment implémenter un MLM dans R ?

**Étape 1.** Construction du modèle *a priori* et exploration des données

###### **Étape 2.** Coder les modèles potentiels et sélection du meilleur modèle

**Étape 3.** Validation du modèle

**Étape 4.** Interprétation et visualisation des résultats

---
# Étape 2. coder le modèle

**Traduisons notre modèle...**

`$$PT_{ijk} \sim Longueur_i + Lac_j + Espèce_k + \varepsilon_{ijk}$$`

**... en code R**

```r
library(lme4)
lmer(Z_TP ~ Z_Length + (1 | Lake) + (1 | Fish_Species),
     data = fish.data, REML = TRUE)
```

* `lmer` ![:faic](arrow-right) fonction "linear mixed model" du package `lme4`
* `(1 | Lake)` ![:faic](arrow-right) indique que les ordonnées à l'origine peuvent varier entre les lacs
* `REML = TRUE` ![:faic](arrow-right) méthode d'estimation

---
# Étape 2. coder le modèle

**Comment faire si on souhaite que les pentes puissent varier?**

---
# Étape 2. coder le modèle

**Différentes structures pour le modèle :**

- `(1 | Lake)` effet aléatoire par lac à l'ordonnée à l'origine
- `(1 + Z_Length | Lake)` effet aléatoire par lac à l'ordonnée à l'origine et la pente en réponse à la longueur corporelle (NB: `(Z_Length | Lake)` donne le même résultat)
- `(-1 + Z_Length | Lake)` pour avoir uniquement l'effet aléatoire sur la pente
- `(1 | Lake) + (1 | Species)`  pour des effets aléatoires croisés (ex: l'effet de l'espèce i ne varie pas entre les lacs)
- `(1 | Lake:Fish_Species)` pour utiliser l'interaction entre 2 facteurs groupant (ex: l'effet de l'espèce i varie d'un lac à l'autre)
- Si votre jeu de données inclus des effets nichés, vous pouvez utiliser `/` pour les déclarer, e.g. `(1 | facteur1 / facteur2)` si `facteur2` est niché dans `facteur1` ([voir ![:faic](stack-exchange)](https://stats.stackexchange.com/questions/228800/crossed-vs-nested-random-effects-how-do-they-differ-and-how-are-they-specified)).

???
Note: S'assurer de bien faire la distinction entre effets croisés (sans interaction), avec interaction et nichés.

---
# Note sur les méthodes d'estimation

REML (*Restricted Maximum Likelihood*) est la méthode par défaut dans la fonction `lmer` (voir `?lmer`).

La méthode de maximum de vraisemblance (ML, pour *Maximum Likelihood*) sous-estime les variances du modèle par un facteur `$(n-k) / n$`, où `$k$` est le nombre d'effets fixes.

La méthode REML corrige ce biais.

Consultez cet [article](https://towardsdatascience.com/maximum-likelihood-ml-vs-reml-78cf79bef2cf) pour plus d'information sur la différence entre ML et REML.

???

Note: Dire que la méthode d'estimation est l'algorithme utiliser pour estimer les paramètres du modèles (ordonnées à l'origine, pentes). Expliquer qu'ils n'ont pas besoin de comprendre l'algorithme en détail, mais qu'ils peuvent suivre le lien pour plus d'informations. La prochaine diapositive résume ce qu'ils devraient savoir.

---
# Note sur les méthodes d'estimation

**On devrait utiliser :**

* **REML** pour comparer des modèles avec différentes structures **d'effets aléatoires** et la même structure d'effets fixes.

* **ML** pour comparer des modèles avec différentes structures **d'effets fixes** et la même structure d'effets aléatoires.

* **ML** pour comparer des modèles **avec effets aléatoires** à des modèles **sans effets aléatoires**.

???

Note: Cela peut sembler complexe, mais les deux prochains défis devraient aider les participants à comprendre les utilisations des difféerentes méthodes d'estimation.

---
# Défi 4 ![:cube]()

Réécrivez le code suivant de façon à ce que les **pentes** de la relation position trophique en fonction de longueur corporelle **varient par lac et par espèce** :

```r
lmer(Z_TP ~ Z_Length + (1 | Lake) + (1 | Fish_Species),
     data = fish.data, REML = TRUE)
# Linear mixed model fit by REML ['lmerMod']
# Formula: Z_TP ~ Z_Length + (1 | Lake) + (1 | Fish_Species)
#    Data: fish.data
# REML criterion at convergence: 72.4662
# Random effects:
#  Groups       Name        Std.Dev.
#  Lake         (Intercept) 0.4516  
#  Fish_Species (Intercept) 0.9301  
#  Residual                 0.2605  
# Number of obs: 180, groups:  Lake, 6; Fish_Species, 3
# Fixed Effects:
# (Intercept)     Z_Length  
#   9.752e-14    4.198e-01
```

---
# Solution ![:cube]()

```r
lmer(Z_TP ~ Z_Length + (1 + Z_Length | Lake) + (1 + Z_Length | Fish_Species),
     data = fish.data, REML = TRUE)
# Linear mixed model fit by REML ['lmerMod']
# Formula: 
# Z_TP ~ Z_Length + (1 + Z_Length | Lake) + (1 + Z_Length | Fish_Species)
#    Data: fish.data
# REML criterion at convergence: 20.5786
# Random effects:
#  Groups       Name        Std.Dev. Corr 
#  Lake         (Intercept) 0.45279       
#               Z_Length    0.02378  -0.82
#  Fish_Species (Intercept) 0.93103       
#               Z_Length    0.15728  1.00 
#  Residual                 0.22341       
# Number of obs: 180, groups:  Lake, 6; Fish_Species, 3
# Fixed Effects:
# (Intercept)     Z_Length  
#  -0.0009025    0.4223738  
# optimizer (nloptwrap) convergence code: 0 (OK) ; 0 optimizer warnings; 1 lme4 warnings
```

---
# Étape 2. sélectionner le meilleur modèle

* Pour déterminer si vous avez construit le meilleur modèle mixte basé sur vos connaissances *a priori*, vous devez comparer ce modèle aux autres modèles alternatifs.

* Avec le jeu de données sur lequel vous travaillez, il y a plusieurs modèles alternatifs qui pourraient mieux correspondre à vos données.

---
# Défi 5 ![:cube]()

Faites une liste de 7 modèles alternatifs qui pourraient être comparés à celui-ci :

```r
lmer(Z_TP ~ Z_Length + (1 | Lake) + (1 | Fish_Species),
     data = fish.data, REML = TRUE)
```

.comment[Note : Si nous avions différents effets fixes entre les modèles ou un modèle sans effets aléatoires, nous aurions eu besoin d'indiquer `REML = FALSE` pour les comparer avec une méthode de vraisemblance comme l'AIC.]

???

Note: Pour savoir si le modèle que vous avez construis basé sur vos connaissances a priori est le meilleur, vous devez le comparer à des modèles avec une structure alternative. Il existe de nombreux modèles alternatifs qui peuvent être construit à partir de la base de données avec laquelle vous travailler.

---
# Solution ![:cube]()

Nous allons aussi construire le **modèle linéaire de base** `lm()` parce qu'il est toujours utile de voir la variation dans les valeurs d'AIC.

```r
M0 <- lm(Z_TP ~ Z_Length, data = fish.data)
```

Par contre, pour comparer ce modèle aux MLMs, il est important de .alert[changer la méthode d'estimation à ML (`REML=FALSE`)] parce que `lm()` n'utilise pas la même méthode d'estimation que `lmer()`.

---
# Solution ![:cube]()

```r
# Modèle linéaire de base
M0 <- lm(Z_TP ~ Z_Length, data = fish.data)
# Modèle complet avec variation des intercepts
M1 <- lmer(Z_TP ~ Z_Length + (1 | Fish_Species) + (1 | Lake), 
 data = fish.data, REML = FALSE)
# Modèle complet avec variation des intercepts et des pentes
M2 <- lmer(Z_TP ~ Z_Length + (1 + Z_Length | Fish_Species) + (1 + Z_Length | Lake),
 data = fish.data, REML = FALSE)
# Pas d'effet lac, les intercepts varient par espèce
M3 <- lmer(Z_TP ~ Z_Length + (1 | Fish_Species), data = fish.data, REML = FALSE)
# Pas d'effet espèces, les intercepts varient par lac
M4 <- lmer(Z_TP ~ Z_Length + (1 | Lake), data = fish.data, REML = FALSE)
# Pas d'effet lac, les intercepts et les pentes varient par espèce
M5 <- lmer(Z_TP ~ Z_Length + (1 + Z_Length | Fish_Species), 
 data = fish.data, REML = FALSE)
# Pas d'effet espèces, les intercepts et les pentes varient par lac
M6 <- lmer(Z_TP ~ Z_Length + (1 + Z_Length | Lake), data = fish.data, REML = FALSE)
# Modèle complet avec variation des intercepts et des pentes par lac
M7 <- lmer(Z_TP ~ Z_Length + (1 | Fish_Species) + (1 + Z_Length | Lake),
 data = fish.data, REML = FALSE)
# Modèle complet avec variation des intercepts et des pentes par espèce
M8 <- lmer(Z_TP ~ Z_Length + (1 + Z_Length | Fish_Species) + (1 | Lake),
 data = fish.data, REML = FALSE)
```

---
# Solution ![:cube]()

Lorsqu'on ajuste des MLM avec `lmer()`, il est possible de faire face à certaines erreurs ou avertissements comme :

* `boundary (singular) fit: see ?isSingular`, voir [cette discussion ![:faic](stack-exchange)](https://stats.stackexchange.com/questions/378939/dealing-with-singular-fit-in-mixed-models)

* `Model failed to converge with max|grad| ...`, voir [cette discussion ![:faic](stack-exchange)](https://stats.stackexchange.com/questions/242109/model-failed-to-converge-warning-in-lmer)

Voici une [liste](https://rdrr.io/cran/lme4/man/troubleshooting.html) de problèmes possibles et la façon de les résoudre.

---
# Étape 2. sélectionner le meilleur modèle

* Maintenant que nous avons une liste de modèles potentiels, nous voulons les comparer entre eux pour sélectionner celui(ceux) qui a(ont) le plus de pouvoir prédictif.

* Les modèles peuvent être comparés en utilisant la fonction `AICc` provenant du package `MuMIn`.

* Le critère d'information Akaike (AIC) est une **mesure de qualité du modèle** pouvant être utilisée pour comparer les modèles.

* `AICc` corrige le biais créé par les faibles tailles d'échantillon.

Pour trouver la valeur AICc d'un modèle, utilisez :

```r
library(MuMIn)
MuMIn::AICc(M1)
# [1] 77.30499
```

---
# Étape 2. sélectionner le meilleur modèle

Pour regrouper toutes les valeurs d'AICc dans un seul tableau, utilisez :

```r
AIC.table <- MuMIn::model.sel(M0, M1, M2, M3, M4, M5, M6, M7, M8)
(AIC.table <- AIC.table[ , c("df", "logLik", "AICc", "delta")])
# df logLik AICc delta
# M8 7 -8.597929 31.84702 0.000000
# M2 9 -8.216019 35.49086 3.643839
# M1 5 -33.480080 77.30499 45.457965
# M7 7 -33.186374 81.02391 49.176890
# M5 6 -128.310995 269.10754 237.260517
# M3 4 -134.532965 277.29450 245.447480
# M4 4 -224.715763 457.66010 425.813076
# M6 6 -224.671201 461.82795 429.980930
# M0 3 -236.861362 479.85909 448.012065
```

* `df` est le degré de liberté
* `logLik` est le log de la vraisemblance
* `delta` est la différence d'AICc avec la valeur la plus petite

Nous avons seulement affiché une partie des résultats retournés par la fonction `model.sel()`, voir `?model.sel` pour plus d'information.

---
# Étape 2. sélectionner le meilleur modèle

Que signifient ces valeurs d'AICc ?

```r
AIC.table
#    df      logLik      AICc      delta
# M8  7   -8.597929  31.84702   0.000000
# M2  9   -8.216019  35.49086   3.643839
# M1  5  -33.480080  77.30499  45.457965
# M7  7  -33.186374  81.02391  49.176890
# M5  6 -128.310995 269.10754 237.260517
# M3  4 -134.532965 277.29450 245.447480
# M4  4 -224.715763 457.66010 425.813076
# M6  6 -224.671201 461.82795 429.980930
# M0  3 -236.861362 479.85909 448.012065
```

* Le modèle avec le plus petit AICc a le plus grand pouvoir prédictif.

* Souvent on considère que deux modèles à +/- 2 unités d'AICc de différence ont un pouvoir prédictif équivalent.

* Examinons de plus proche M8 and M2. On peut exclure les autres modèles, car ils ont des AICc beaucoup plus élevés.

---
# Étape 2. sélectionner le meilleur modèle

```r
M8 <- lmer(Z_TP ~ Z_Length + (1 + Z_Length | Fish_Species) + (1 | Lake),
 data = fish.data, REML = TRUE)

M2 <- lmer(Z_TP ~ Z_Length + (1 + Z_Length | Fish_Species) + (1 + Z_Length | Lake),
 data = fish.data, REML = TRUE)

MuMIn::model.sel(M2,M8)[ , c("df", "logLik", "AICc", "delta")]
#    df    logLik     AICc    delta
# M8  7 -10.84011 36.33137 0.000000
# M2  9 -10.28932 39.63747 3.306098
```

Le modèle `M8` semble être le meilleur modèle parmi ceux qu'on a testés.

Notez qu'on utilise maintenant REML (i.e. `REML=TRUE`) étant donné qu'on compare deux modèles avec des effets aléatoires nichés et avec la même structure d'effets fixes.

---
# Étape 2. sélectionner le meilleur modèle

Quelle est la structure du meilleur modèle?

```r
M8 <- lmer(Z_TP ~ Z_Length + (1 + Z_Length | Fish_Species) + (1 | Lake),
 data = fish.data, REML = TRUE)
```

Les ordonnées à l'origine et les pentes de l'effet de la longueur sur la position trophique peuvent varier selon l'espèce de poissons, mais seulement des ordonnées à l'origine peuvent varier par lac.

---
# Étape 2. sélectionner le meilleur modèle

Une fois que le meilleur modèle est sélectionné, il faut remettre la méthode d'estimation à `REML = TRUE`.

```r
M8 <- lmer(Z_TP ~ Z_Length + (1 + Z_Length | Fish_Species) + (1 | Lake),
 data = fish.data, REML = TRUE)
```

---
exclude: true

# Défi 5 ![:cube]()

Prenez 2 minutes avec votre voisin pour étudier la structure du modèle M2.

Comment diffère-t-il de M8 d'un point de vue écologique?

Pourquoi n'est-il pas surprenant que sa valeur d'AICc était la deuxième meilleure?

```r
M8 <- lmer(Z_TP ~ Z_Length + (1 + Z_Length | Fish_Species) + (1 | Lake),
 data = fish.data, REML = TRUE)

M2 <- lmer(Z_TP ~ Z_Length + (1 + Z_Length | Fish_Species) + (1 + Z_Length | Lake),
 data = fish.data, REML = TRUE)
```

---
exclude: true

# Solution ![:cube]()

**Discussion de groupe...**

--
exclude: true

.alert[M2] La position trophique est en fonction de la longueur. L'intercept et l'effet de la longueur sur la position trophique peuvent varier selon l'espèce de poissons et le lac.

* .small[les facteurs intrinsèques des espèces et des lacs sont à la base des relations différentes entre la position trophique et la longueur (i.e. pentes et intercepts)]

.alert[M8] La position trophique est en fonction de la longueur. L'intercept et l'effet de la longueur sur la position trophique peut varier selon l’espèce de poissons, mais seulement l'intercept peut varier par lac.

* .small[seulement les facteurs intrinsèques des espèces sont responsables des différentes relations (i.e. pentes) et en moyenne, les positions trophiques pourraient être supérieures ou inférieures d’un lac à l’autre (e.g. intercepts).]

---
# Comment implémenter un MLM dans R ?

**Étape 1.** Construction du modèle *a priori* et exploration des données

**Étape 2.** Coder les modèles potentiels et sélection du meilleur modèle

###### **Étape 3.** Validation du modèle

**Étape 4.** Interprétation et visualisation des résultats

---
# Étape 3. validation du modèle

**Vous devez vérifier que le modèle respecte toutes les suppositions de base :**

**3.1 Vérifier l'homogénéité de la variance**
  - Graphique des valeurs prédites en fonction des valeurs résiduelles

**3.2 Vérifier l'indépendance des résidus**
  - Graphique des résidus vs chaque covariable du modèle
  - Graphique des résidus vs chaque covariable non incluse dans le modèle

**3.3 Vérifier la normalité**
  - Histogramme des résidus

---
# Étape 3. validation du modèle

**3.1 Vérifier l'homogénéité de la variance**

```r
plot(resid(M8) ~ fitted(M8), 
     xlab = 'Valeurs prédites', 
     ylab = 'Résidus normalisés')
abline(h = 0, lty = 2)
```

Étendue homogène des résidus ![:faic](arrow-right) la supposition est respectée!

---
# Étape 3. validation du modèle

**3.1 Vérifier l'homogénéité de la variance**

???

Note: L'hétéroscédasticité est lorsque la variance des résidus du modèle (i.e., erreur) varie avec une variable explicative. Lorsque l'on graphe les résidus avec une variable explicative, nous observons une forme conique des résidus comme les examples ci-dessous. L'inverse d'hétéoscedasticité est l'homoscédasticité. L'homoscédasticité est lorsque la variance des résidus ne varie pas avec les variables explicatives. C'est ce qu'on veut!
---
# Étape 3. validation du modèle

**3.2 Vérifier l'indépendance des résidus avec chaque covariable**

```r
plot(resid(M8) ~ fish.data$Z_Length, 
     xlab = "Longueur", ylab = "Résidus normalisés")
abline(h = 0, lty = 2)

boxplot(resid(M8) ~ Fish_Species, data = fish.data, 
        xlab = "Espèces", ylab = "Résidus normalisés")
abline(h = 0, lty = 2)

boxplot(resid(M8) ~ Lake, data = fish.data, 
        xlab = "Lacs", ylab = "Résidus normalisés")
abline(h = 0, lty = 2)
```

---
# Étape 3. validation du modèle

**3.2 Vérifier l'indépendance des résidus avec chaque covariable**

Étendue homogène des résidus autour de 0 ![:faic](arrow-right) pas de patron des résidus en fonction de la variable, la supposition est respectée!

.comment[Note : Les regroupements de données sont causés par la structure des données, où des poissons de seulement 5 classes de taille  (grand, petit, et trois groupes entre les deux) étaient capturés.]

---
# Étape 3. validation du modèle

**3.2 Vérifier l'indépendance des résidus avec chaque covariable**

- Graphique des résidus vs chaque covariable non incluse dans le modèle

- Si vous observez des patrons dans ce graphique, vous saurez qu'il y a de la variation dans votre jeu de données qui pourrait être expliquée par ces covariables. Vous devriez considérer d'inclure ces variables dans votre modèle.

- Étant donnée que nous avons inclus toutes les variables mesurées dans notre modèle, nous ne pouvons pas faire cette étape.

---
# Étape 3. validation du modèle

**3.3 Vérifier la normalité des résidus**

* Des résidus suivant une distribution normale indiquent que le modèle n'est pas biaisé.

```r
hist(resid(M8))
```

---
# Comment implémenter un MLM dans R ?

**Étape 1.** Construction du modèle *a priori* et exploration des données

**Étape 2.** Coder les modèles potentiels et sélection du meilleur modèle

**Étape 3.** Validation du modèle

###### **Étape 4.** Interprétation et visualisation des résultats

---
# Étape 4. interprétation et visualisation

```r
(summ_M8 <- summary(M8))
# Linear mixed model fit by REML ['lmerMod']
# Formula: Z_TP ~ Z_Length + (1 + Z_Length | Fish_Species) + (1 | Lake)
# Data: fish.data
# 
# REML criterion at convergence: 21.7
# 
# Scaled residuals: 
# Min 1Q Median 3Q Max 
# -2.77187 -0.60166 0.05589 0.64239 2.27776 
# 
# Random effects:
# Groups Name Variance Std.Dev. Corr
# Lake (Intercept) 0.20504 0.4528 
# Fish_Species (Intercept) 0.86715 0.9312 
# Z_Length 0.02466 0.1570 1.00
# Residual 0.05039 0.2245 
# Number of obs: 180, groups: Lake, 6; Fish_Species, 3
# 
# Fixed effects:
# Estimate Std. Error t value
# (Intercept) -0.0009059 0.5687733 -0.002
# Z_Length 0.4222697 0.0922117 4.579
# 
# Correlation of Fixed Effects:
# (Intr)
# Z_Length 0.929 
# optimizer (nloptwrap) convergence code: 0 (OK)
# boundary (singular) fit: see help('isSingular')
```

???

Note: Prenons un moment pour explorer un peu plus notre modèle final à l'aide de la fonction summary(). Comment pouvons-nous interpréter ces informations?

---
# Étape 4. interprétation et visualisation

# Random effects:
    #  Groups       Name        Variance Std.Dev. Corr
    #  Lake         (Intercept) 0.20500  0.4528
    #  Fish_Species (Intercept) 0.86621  0.9307
    #               Z_Length    0.02464  0.1570   1.00
    #  Residual                 0.05040  0.2245

- `Groups`: facteurs de regroupement
- `Name`:
  - `(Intercept)` pour l'ordonnée à l'origine,
  - ou le nom de la variable sur lequel porte l'effet mixte dans le cas d'une pente aléatoire (`Z_length` dans notre exemple)
- `Variance`: variance estimée de l'effet (`Std.Dev.` est l'écart type de cette valeur)
- `Corr`: corrélation entre la pente aléatoire et l'ordonnée à l'origine aléatoire pour un groupement donné (voir [cette discussion ![:faic](stack-exchange)](https://stats.stackexchange.com/questions/320978/understanding-and-coding-random-intercept-correlation-lmer))

---
# Étape 4. interprétation et visualisation

# Fixed effects:
    #              Estimate Std. Error t value
    # (Intercept) -0.000906   0.568493  -0.002
    # Z_Length     0.422270   0.092170   4.581

Cette partie présente l'estimation des effets fixes. Une valeur de la statistique T [(test de Student)](https://en.wikipedia.org/wiki/T-statistic) est retournée **sans p-value** (c'est un choix des auteurs du package, voir pourquoi dans [cette discussion](https://stats.stackexchange.com/questions/185360/t-value-associated-with-nlme-lme4)).

Cette statistique peut être utilisée telle quelle. Vous pouvez aussi calculer
l’intervalle de confiance (IC) à 95% avec cette table en utilisant :

$$ IC = Estimate \pm 1.96*Std.Error $$
Si 0 est dans cet intervalle, alors le paramètre n’est pas significativement
différent de zéro au seuil `$\alpha$` = 0.05.

???

Note: Passer section par section pour bien comprendre ce qui est présenté.

La sortie est divisée en deux: la description des facteurs aléatoires et les facteurs fixes.

---
# Étape 4. interprétation et visualisation

**Quelques fonctions utiles**

- `coef(M8)` et `ranef(M8)` retournent les effets aléatoires du modèle M8

- `coef(summary(M8))` retourne les effets fixes

- `sigma(M8)` retourne l’écart type des résidus

- `fitted(M8)` retourne les valeurs prédites par le modèle

- `residuals(M8)` retourne les résidus

---
# Défi 6 ![:cube]()

1. Quelle est la pente et son intervalle de confiance de la variable Z_Length dans le modèle M8?

2. Est-ce que la pente de Z_Length est significativement différente de 0 ?

---
# Solution ![:cube]()

1. Quelle est la pente et son intervalle de confiance de la variable Z_Length dans le modèle M8?

- pente = 0.422;

- limite supérieure de l’IC = 0.4223 + 0.09*1.96 = 0.5987

- limite inférieure de l’IC = 0.4223 - 0.09*1.96 = 0.2459

2. Est-ce que la pente de Z_Length est significativement différente de 0 ?

- Oui, car l'IC [0.2459, 0.5987] n'inclut pas 0

???

Note: L'intervalle de confiance de la pente n'inclue pas 0 (0 = aucun effet), la pente est donc significativement différente de 0. Cela signifie que la longueur du poisson influence sa position trophique (un plus long poisson a, en moyenne, un plus haut niveau trophique). Est-ce que ça fait du sens? Oui, les gros poissons mangent les petits poissons.

---
# Défi 7 ![:cube]()

Il est possible de visualiser graphiquement les différentes intercepts et pentes du modèle pour mieux interpréter les résultats

Prenez 2 minutes pour réfléchir aux différentes façons de représenter les résultats de M8.

*Indice: considérez les différents "niveaux" du modèle*

---
# Solution ![:cube]()

a) Figure avec toutes les données regroupées

b) Figure par espèce

c) Figure par lac

---
# Solution ![:cube]()

Pour faire ces figures, il nous faut:

- Les coefficients du modèle complet qui sont dans le résumé du modèle

```r
summ_M8$coefficients
#                  Estimate Std. Error      t value
# (Intercept) -0.0009058974 0.56877327 -0.001592722
# Z_Length     0.4222697238 0.09221166  4.579352788
```

- Intercept = `$-9.0589745\times 10^{-4}$`
- Pente = `$0.4222697$`

---
# Solution ![:cube]()

Pour faire ces figures, il nous faut:

- Les coefficients pour chaque niveau du modèle qu'on obtient avec la fonction `coef`

```r
coef(M8)
# $Lake
#     (Intercept)  Z_Length
# L1 -0.085984071 0.4222697
# L2  0.002205209 0.4222697
# L3 -0.301816557 0.4222697
# L4 -0.574039728 0.4222697
# L5  0.218650140 0.4222697
# L6  0.735549622 0.4222697
# 
# $Fish_Species
#    (Intercept)  Z_Length
# S1  -1.0752985 0.2410746
# S2   0.5597871 0.5168300
# S3   0.5127938 0.5089046
# 
# attr(,"class")
# [1] "coef.mer"
```

---
# Solution ![:cube]()

a) Figure avec toutes les données regroupées

```r
library(ggplot2)

# Thème ggplot simplifié
fig <- theme_bw() +
 theme(panel.grid.minor=element_blank(), 
 panel.grid.major=element_blank(),
 panel.background=element_blank()) +
 theme(strip.background=element_blank(), 
 strip.text.y = element_text()) +
 theme(legend.background=element_blank()) +
 theme(legend.key=element_blank()) +
 theme(panel.border = element_rect(colour = "black", fill=NA))

plot <- ggplot(aes(Z_Length, Z_TP), data = fish.data)
Plot_AllData <- plot + geom_point() +
 xlab("Longueur (mm)") + 
 ylab("Position trophique") +
 labs(title = "Toutes les données") + fig

Plot_AllData + geom_abline(intercept = summ_M8$coefficients[1,1], 
                           slope     = summ_M8$coefficients[2,1])
```

---
# Solution ![:cube]()

a) Figure avec toutes les données regroupées
<img src="workshop07-pres-fr_files/figure-html/unnamed-chunk-39-1.png" width="432" style="display: block; margin: auto;" />

---
# Solution ![:cube]()

b) Figure par espèce

```r
# mettre les coefs dans un tableau pour les rendre plus faciles à manipuler
Lake.coef <- coef(M8)$Lake
colnames(Lake.coef) <- c("Intercept", "Slope")
Species.coef <- coef(M8)$Fish_Species
colnames(Species.coef) <- c("Intercept", "Slope")

Plot_BySpecies <- plot + 
 geom_point(aes(colour = factor(Fish_Species)), size = 4) +
 xlab("Longueur (mm)") + ylab("Position trophique") +
 labs(title = "Par espèce") + fig

# Ajoutez les lignes de régression pour chaque espèce
Plot_BySpecies +
  geom_abline(intercept = Species.coef[1,1], 
              slope     = Species.coef[1,2], col = "coral2") +
  geom_abline(intercept = Species.coef[2,1], 
              slope     = Species.coef[2,2], col = "green4") +
  geom_abline(intercept = Species.coef[3,1], 
              slope     = Species.coef[3,2], col = "blue1")

```

---
# Solution ![:cube]()

b) Figure par espèce

---
# Solution ![:cube]()

c) Figure par lac

```r

Plot_ByLake <- plot + 
 geom_point(aes(colour = factor(Lake)), size = 4) +
 xlab("Length (mm)") + ylab("Trophic Position") +
 labs(title = "By Lake") + fig

# Ajouter les lignes de régression avec les intercepts spécifiques à chaque lac
Plot_ByLake +
  geom_abline(intercept = Lake.coef[1,1], 
              slope     = Lake.coef[1,2], col = "coral2") +
  geom_abline(intercept = Lake.coef[2,1], 
              slope     = Lake.coef[2,2], col = "khaki4") +
  geom_abline(intercept = Lake.coef[3,1], 
              slope     = Lake.coef[3,2], col = "green4") +
  geom_abline(intercept = Lake.coef[4,1], 
              slope     = Lake.coef[4,2], col = "darkgoldenrod") +
  geom_abline(intercept = Lake.coef[5,1], 
              slope     = Lake.coef[5,2], col = "royalblue1") +
  geom_abline(intercept = Lake.coef[6,1], 
              slope     = Lake.coef[6,2], col = "magenta3")

```

---
# Solution ![:cube]()

c) Figure par lac
<img src="workshop07-pres-fr_files/figure-html/unnamed-chunk-43-1.png" width="576" style="display: block; margin: auto;" />

---

# Modèle mixtes et données en écologie

Les modèles mixtes sont très utiles pour prendre en compte la structure complexe des données en écologie tout en permettant de ne pas perdre beaucoup de degrés de liberté

---
exclude: true
# Défi 8 ![:cube]()

**Situation:**

* Vous avez inventorié la richesse **dans 1000 quadrats** qui sont dans **10 sites différents** qui sont également dans **10 forêts différentes**.

* Vous avez de plus **mesuré la productivité** dans chaque **quadrat**.

* Vous désirez savoir si la productivité est un bon prédicteur de biodiversité.

---
exclude: true
# Solution ![:cube]()

```r
lmer(Biodiv ~ Productivite + (1 | Foret / Site))
```

Ici les effets aléatoires sont nichés (i.e. Sites dans forêt) et non croisés.

Pourquoi utiliser `(1 | Foret / Site)` plutôt que `(1 | Foret) + (1 | Site)` ? Regardez [cette réponse sur ![:faic](stack-exchange)](https://stats.stackexchange.com/questions/228800/crossed-vs-nested-random-effects-how-do-they-differ-and-how-are-they-specified)!

---
exclude: true

# Défi 8 ![:cube]()

**Situation:**

* Vous avez récolté **200 poissons** dans **12 sites différents** distribués également dans **4 habitats** différents qui se retrouvent dans **un même lac**.

* Vous avez mesuré la **longueur de chaque poisson** et la **quantité de mercure dans ses tissus**.

* Vous désirez savoir si l'habitat est un bon prédicteur de la concentration en mercure.

---
exclude: true

# Solution ![:cube]()

```r
lmer(Mercure ~ Longueur * Habitat + (1 | Site))
```

---
exclude: true

# Défi 9 ![:cube]()

* Discutez du jeu de données sur lequel vous travaillez avec votre voisin et déterminez si un modèle mixte serait approprié.

* Si oui, travaillez ensemble pour écrire le code que vous utiliseriez pour faire ce modèle dans R.

* Si non, imaginez un jeu de données fictif pour lequel un modèle mixte serait approprié et codez ce modèle.

---
class: inverse, center, middle

# GLMMs

---
# Modèles Linéaires Généralisés Mixtes (GLMMs)

**Extension des GLMs tenant compte de structures supplémentaires dans les données**

Ils suivent les étapes similaires à celles introduites lors de la section sur les LMMs :

**1.** Incorporent les effets aléatoires (comme les LMMs)

**2.** Permettent de gérer des données non-normales (en laissant les erreurs prendre différentes familles de distribution - e.g Poisson ou binomial négatif) (comme les GLMs)

---
# Comment modéliser un GLMM avec R

Chargez les données `Arabidopsis` `arabidopsis.csv` dans `R`.

Cette base de données décrit la réponse génétique de l'espèce `Arabidopsis` à la fertilisation et l'herbivorie.

Pour plus de détails, vous pouvez consulter [ceci](https://www.rdocumentation.org/packages/lme4/versions/1.1-27.1/topics/Arabidopsis).

```r
# Voici la description des colonnes de ce jeu de données:

# reg: facteur avec un niveau pour chaque région
# popu: facteur avec un niveau pour chaque population
# gen: facteur avec un niveau pour chaque génotype
# nutrient: facteur avec niveau bas (valeur = 1) ou haut (valeur = 8)
# amd: facteur précisant l'absence ou la présence d'herbivorie
# total.fruits: nombre entier indiquant le nombre de fruits par plante
```

L'effet de la disponibilité de nutriments et d'herbivorie (**effets fixes**) sur la production de fruits d'*Arabidopsis thaliana* (**variable réponse**) a été évalué en mesurant 625 plantes à travers neuf populations différentes, constituées chacune de 2 à 3 génotypes (**effets aléatoires**).

---
# Choisir la distribution des erreurs

La variable réponse constitue des données d'abondance, donc nous devons choisir une **distribution de Poisson** (i.e variance égale à la moyenne).

Vous pouvez lire plus sur la distribution des erreurs [ici](http://ganymede.nmsu.edu/holtz/a535/ay535notes/node3.html).

Cependant, comme nous le verrons, la variance de chaque groupe augmente beaucoup plus rapidement que prévu... Cela indique que les données sont surdispersées.

---
# Surdispersion: rappel

Pour une distribution de Poisson, `$var[y] = μ$`. Cependant, en pratique, on constate que la variance des données dépasse souvent `$μ$`, indiquant une *surdispersion* dans les paramètres du modèle.

La surdispersion est due au fait que la moyenne `$μ$` varie intrinsèquement,
même lorsque toutes les variables explicatives sont fixes, ou parce que les événements qui sont comptés sont corrélés positivement.

**Pourquoi devrions-nous nous intéresser à la surdispersion?**

> Les tests sur les variables explicatives apparaîtront généralement plus significatifs et les intervalles de confiance des paramètres seront plus étroits que les données ne le justifieraient!

---
# Exploration de la variance

Pour illustrer l'hétérogénéité de la variance, nous allons créer des boîtes à moustaches (boxplots) du **log** du nombre total de fruits (**variable réponse**) par rapport aux différents facteurs environnementaux.

Créons d'abord de nouvelles variables qui représentent toutes les combinaisons de **nutriments** x **herbivorie** x **facteur aléatoire**.

```r
dat.tf <- within(dat.tf,
{
 # génotype x nutriments x herbivorie
 gna <- interaction(gen, nutrient, amd)
 gna <- reorder(gna, total.fruits, mean)
 # population x nutriments x herbivorie
 pna <- interaction(popu, nutrient, amd)
 pna <- reorder(pna, total.fruits, mean)
})
```

---
# Exploration de la variance

```r
# Boxplot du total des fruits vs interaction génotype x nutriments x herbivorie
ggplot(data = dat.tf, aes(factor(x = gna), y = log(total.fruits + 1))) +
  geom_boxplot(colour = "skyblue2", outlier.shape = 21,
  outlier.colour = "skyblue2") +
  ylab("log (Fruits totaux)\n") + # \n créé un espace après le titre
  xlab("\nGénotype x nutriments x herbivorie") + # espace avant le titre
  theme_bw() + theme(axis.text.x = element_blank()) +
  stat_summary(fun = mean, geom = "point", colour = "red")
```

<img src="workshop07-pres-fr_files/figure-html/unnamed-chunk-50-1.png" width="576" style="display: block; margin: auto;" />
]

.comment[De même, la variance du total de fruits montre une grande hétérogénéité entre les populations (population x nutriments x herbivorie).]

---
# Choisir la distribution des erreurs

Comme nous venons de le voir, il existe une importante hétérogénéité parmi la variance de chaque groupe, et ce, même lorsque la variable réponse est transformée (i.e. log).

Si nous représentons graphiquement les **écarts vs moyennes par groupe** (génotype x nutriments x herbivorie), on voit que la distribution de Poisson est la moins appropriée (car les écarts augmentent beaucoup plus rapidement que la moyenne).

---
# GLMM Poisson

**Compte tenu de la relation moyenne-variance, nous avons besoin d'un modèle qui tient compte de la surdispersion.**

Pour comprendre pourquoi, commençons avec un modèle avec une distribution de Poisson.

Pour lancer un GLMM dans `R`, nous faisons appel à la fonction `glmer()` de la librairie `lme4` :

```r
library(lme4)
mp1 <- glmer(total.fruits ~ nutrient*amd + rack + status +
 (1|popu)+
 (1|gen),
 data = dat.tf, family = "poisson")
```

**Effets aléatoires** : `(1|popu)` et `(1|gen)`. Nous faisons varier les ordonnées à l'origine pour les différentes populations (`popu`) et génotypes (`gen`).

---
# Vérification de la surdispersion

On vérifie la surdispersion en utilisant la fonction `overdisp_fun()` (Bolker *et al*. 2011) qui divise la déviance des résidus (de Pearson) par leurs degrés de liberté.

La fonction évalue si le **rapport est plus grand que 1**.

Effectuons ce test :

```r
# Téléchargez le code glmm_funs.R de la page wiki et sourcez le pour exécuter la fonction dans R
source(file = "data/glmm_funs.R")
# Surdispersion?
overdisp_fun(mp1)
# chisq ratio p logp 
# 15755.86835 25.57771 0.00000 -6578.47028
```
--
.alert[Le rapport est significativement > 1]

Comme on s'y attendait, nous devons modéliser une **distribution différente** où la variance augmente plus rapidement que la moyenne.

---
# GLMM binomiale negative .small[(Poisson-gamma)]

**Option 1.** La distribution **binomiale négative** satisfait la supposition que la variance est proportionnelle au carré de la moyenne.

--
On peut modéliser cette distribution avec la fonction `glmer.nb()` :

```r
mnb1 <- glmer.nb(total.fruits ~ nutrient*amd + rack + status +
 (1|popu) +
 (1|gen),
 data = dat.tf, control = glmerControl(optimizer = "bobyqa"))
# Control spécifie la façon dont nous optimisons les valeurs des paramètres
```
--
On évalue à nouveau la présence de surdispersion :

```r
# Surdispersion?
overdisp_fun(mnb1)
#         chisq         ratio             p          logp 
# 721.034461131   1.170510489   0.002143425  -6.145350288
```

---
# GLMM Poisson-lognormale

**Option 2.** La distribution **Poisson-lognormale**

On réalise cette distribution en ajoutant un **effet aléatoire de niveau d'observation** dans le modèle.

Voir Harrison (2014) pour plus de détails https://doi.org/10.7717/peerj.616.

--
Pour ce faire, nous créons tout d'abord une variable nommée `X` :

```r
# Cette variable est déjà dans vos données "dat.tf", mais voici comment la créer.
dat.tf$X <- 1:nrow(dat.tf)
```

On traite la surdispersion en ajoutant l'effet `(1|X)` dans la formule :

```r
mpl1 <- glmer(total.fruits ~ nutrient*amd + rack + status +
 (1|X) +
 (1|popu) +
 (1|gen),
data = dat.tf, family = "poisson",
control = glmerControl(optimizer = "bobyqa"))
```

---
# GLMM Poisson-lognormale

On évalue finalement la présence de surdispersion :

```r
overdisp_fun(mpl1)
#         chisq         ratio             p          logp 
#  1.775362e+02  2.886767e-01  1.000000e+00 -3.755645e-73
```

---
# GLMM Poisson-lognormale

**Représentons graphiquement les paramètres du modèle**

Peut être obtenu en utilisant la fonction `coefplot2()`de la librairie `coefplot2` :

.alert[![:faic](warning) Cette librairie n'est pas sur le CRAN! On utilise la librairie `remotes` pour l'installer à partir de GitHub]

```r
if (!require("coefplot2"))
  remotes::install_github("palday/coefplot2", 
                          subdir = "pkg", 
                          upgrade = "always", 
                          quiet = TRUE)

library(coefplot2)
```

---
# GLMM Poisson-lognormale

```r
# Effets aléatoires
coefplot2(mpl1, ptype = "vcov", intercept = TRUE, main = "Variance des effets aléatoires")
```

<img src="workshop07-pres-fr_files/figure-html/unnamed-chunk-59-1.png" width="432" style="display: block; margin: auto;" />
]
.pull-right[

```r
# Effets fixes
coefplot2(mpl1, intercept = TRUE, main = "Coefficients des effets fixes")
```

<img src="workshop07-pres-fr_files/figure-html/unnamed-chunk-60-1.png" width="432" style="display: block; margin: auto;" />
]

.comment[Note : barres d'erreur visibles seulement pour les effets fixes, car `glmer()` ne modélise pas d'incertitude pour les effets aléatoires.]

---
# GLMM Poisson-lognormale

**Visualisation des effets aléatoires**

Vous pouvez extraire les effets aléatoires en utilisant la fonction `ranef()` et les tracer en utilisant un `dotplot()` de la librairie `lattice`.

--
On constate une variation **inter-population** :

- Les populations espagnoles (SP) ont des valeurs plus élevées que les populations suédoises (SW) et néerlandaises (NL)

On constate une faible variation **inter-génotype** :
- Les différences entre les génotypes semblent induites par le génotype 34

```r
library(gridExtra)
library(lattice)
# dotplot code
pp <- list(layout.widths = list(left.padding = 0, right.padding = 0),
 layout.heights = list(top.padding = 0, bottom.padding = 0))
r2 <- ranef(mpl1, condVar = TRUE)
d2 <- dotplot(r2, par.settings = pp)
grid.arrange(d2$gen, d2$popu, nrow = 1)
```

---
# GLMM Poisson-lognormale

**Visualisation des effets aléatoires**

---
# Sélection du modèle

Les mêmes méthodes peuvent être utilisées avec un glmm ou lmm pour choisir entre des modèles avec différentes ordonnées à l'origine et/ou pentes aléatoires ainsi que pour choisir les effets fixes à conserver dans le modèle final.

En voici deux :

- une **approche de la théorie de l'information** (e.g., AICc - Atelier 5)
- une **approche fréquentiste** (où l'importance de chaque terme est évaluée en faisant un test de rapport de vraisemblance; LRT, avec la fonction `anova()`)

---
# Sélection du modèle

Nous dérivons d'abord les modèles potentiels et les comparons en utilisant l'AICc.comment[*] :

```r
mpl2 <- update(mpl1, . ~ . - rack) # modèle sans rack
mpl3 <- update(mpl1, . ~ . - status) # modèle sans status
mpl4 <- update(mpl1, . ~ . - amd:nutrient) # modèle sans interaction amd:nutrient

aic_tab <- MuMIn::model.sel(mpl1, mpl2, mpl3, mpl4)
(round(aic_table <- aic_tab[ , c("AICc", "delta", "df")], digits = 2))
# AICc delta df
# mpl1 5015.73 0.00 10
# mpl4 5017.11 1.38 9
# mpl3 5017.22 1.49 8
# mpl2 5070.75 55.02 9
```

.comment[*NB: Nous ne couvrons pas tous les modèles possibles ci-dessus, cependant, l'interaction `amd:nutriments` ne peut être évaluée que si les effets additifs de amd et de nutriments sont présents dans le modèle.
]

---
# Sélection du modèle

Nous pouvons aussi utiliser les fonctions `drop1()` et `dfun()` pour évaluer nos effets fixes (`dfun()` convertit les valeurs AIC retournées par `drop1()` en valeurs `$\Delta$`AIC)

```r
dd_LRT <- drop1(mpl1, test = "Chisq")
(dd_AIC <- dfun(drop1(mpl1)))
# Single term deletions
# 
# Model:
# total.fruits ~ nutrient * amd + rack + status + (1 | X) + (1 | 
# popu) + (1 | gen)
# npar dAIC
# <none> 0.000
# rack 1 55.083
# status 2 1.612
# nutrient:amd 1 1.444
```
]

- Fort effet de **rack** (dAIC = 55.08 si on enlève cette variable)
- Effets de **status** et de l'**interaction** sont faibles (dAIC < 2)

---
# Sélection du modèle
Commençons par **enlever l'interaction non significative** afin de tester les effets principaux de nutriments et d'herbivorie.

```r
mpl2 <- update(mpl1, . ~ . - amd:nutrient)

mpl3 <- update(mpl2, . ~ . - rack) # pas de rack ou d'interaction
mpl4 <- update(mpl2, . ~ . - status) # pas de status ou d'interaction
mpl5 <- update(mpl2, . ~ . - nutrient) # pas de nutrient ou d'interaction
mpl6 <- update(mpl2, . ~ . - amd) # pas d'herbivorie ou d'interaction
```
---
# Sélection du modèle
Sélectionnez le meilleur modèle avec la méthode de votre choix :

**Méthode par `AICc`**

```r
aic_tab2 <- MuMIn::model.sel(mpl2, mpl3, mpl4, mpl5, mpl6)
(round(aic_table2 <- aic_tab2[ , c("AICc", "delta", "df")], digits = 2))
# AICc delta df
# mpl2 5017.11 0.00 9
# mpl4 5018.28 1.17 7
# mpl6 5027.27 10.16 8
# mpl3 5071.28 54.17 8
# mpl5 5152.74 135.63 8
```

**Méthode avec `drop1()`**

```r
dd_LRT2 <- drop1(mpl2, test = "Chisq")
dd_AIC2 <- dfun(drop1(mpl2))
```

---
# Sélection du modèle
**Quels sont nos constats ?**

Les **nutriments** et l'**herbivorie** ont un effet important (grand changement d'AIC de `$135.6$` (`mpl5`) et `$10.2$` (`mpl6`) si l'un ou l'autre sont supprimés, respectivement)

**Notre modèle final inclut donc :**

Effets fixes
* nutriments 
* herbivorie
* la variable nuisance de rack.

Effets aléatoires 
* effet du niveau d'observation `(1|X)` 
* populations `(1|popu)` 
* génotypes `(1|gen)`

---
# Défi 8 ![:cube]()

En utilisant l'ensemble de données `inverts` (temps de développement larvaire (`PLD`) de 74 espèces d'invertébrés et vertébrés marins élevés à différentes températures et temps), répondez aux questions suivantes:

- Quel est l'effet du type d'alimentation et du climat (**effets fixes**) sur `PLD`?
- Est-ce que cette relation varie selon les taxons (**effets aléatoires**)?
- Quelle est la **meilleure famille de distributions** pour ces données?
- Finalement, une fois que vous avez déterminé la meilleure famille de distribution, réévaluez vos effets fixes et aléatoires.

---
# Solution ![:cube]()

```r
inverts <- read.csv('data/inverts.csv', header = TRUE)
head(inverts)
table(inverts$temp, inverts$feeding.type)

mod.glm <- glm(PLD ~ temp + feeding.type, family = poisson(), data = inverts)
summary(mod.glm)
drop1(mod.glm, test = "Chisq")

boxplot(PLD ~ temp,  data = inverts)
boxplot(PLD ~ feeding.type ,  data = inverts)

boxplot(predict(mod.glm, type = "response")~inverts$temp)

modglm <- glm(PLD ~ temp + feeding.type, family = poisson(), data = inverts)
```

---
# Solution ![:cube]()

```r
r2 <- ranef(mpl1, condVar = TRUE)
d2 <- dotplot(r2, par.settings = pp)

plot(aggregate(PLD ~ taxon, FUN = mean, data = inverts)[,2], aggregate(PLD ~ taxon, FUN = var, data = inverts)[,2], pch = 19)
abline(a = 0, b = 1, lty = 2)

mod.glmer <- glmer.nb(PLD ~ temp + feeding.type + (1|taxon), data = inverts)
mod.glm <- glm.nb(PLD ~ temp + feeding.type, data = inverts)

plot(aggregate(PLD ~ taxon, FUN = var, data = inverts)[,2], aggregate(PLD ~ taxon, FUN = mean, data = inverts)[,2])
abline(a = 0, b = 1, lty = 2 )
```

---
# Ressources additionnelles

.center[
![:scale 16%](images/book1.jpg) ![:scale 18%](images/book2.jpg) ![:scale 16%](images/book3.jpg)
]

**Libraries connues pour des (G)LMM**

* Fréquentiste : `nlme`, `lme4`, `glmmTMB`
* Bayésien : `brms`, `rstan`, `rstanarm`, `MCMCglmm`

**Articles**
* [Harrison *et al.* (2018), *PeerJ*, DOI 10.7717/peerj.4794
](http://dx.doi.org/10.7717/peerj.4794)
* [Silk *et al.* (2020), *PeerJ*,  DOI 10.7717/peerj.9522
](https://doi.org/10.7717/peerj.9522)
* [Schielzeth *et al.* (2020), *Methods Ecol Evol.*,  DOI: 10.1111/2041-210X.13434
]( https://doi.org/10.1111/2041-210X.13434)
* [Zuur & Ieno (2016), *Methods Ecol Evol.*,  DOI: 10.1111/2041-210X.12577
]( https://doi.org/10.1111/2041-210X.12577)

---
class: inverse, center, bottom

# Merci de votre participation à cet atelier!

![:scale 50%](images/qcbs_logo.png)