Atelier 6: Modèles linéaires généralisés

class: center, middle, inverse, title-slide

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

---

class: inverse, center, middle

# À propos de cet atelier

[![badge](https://img.shields.io/static/v1?style=for-the-badge&label=Présentation&message=06&color=BF616A)](https://r.qcbs.ca/workshop06/pres-fr/workshop06-pres-fr.html)    [![badge](https://img.shields.io/static/v1?style=for-the-badge&label=Script&message=06&color=D08770&logo=r)](https://r.qcbs.ca/workshop06/book-fr/workshop06-script-fr.R) 
[![badge](https://img.shields.io/static/v1?style=for-the-badge&label=Livre&message=06&color=EBCB8B)](https://r.qcbs.ca/workshop06/book-fr/index.html)
[![badge](https://img.shields.io/static/v1?style=for-the-badge&label=Site&message=r.qcbs.ca&color=A3BE8C)](https://r.qcbs.ca/fr/workshops/r-workshop-06/)
[![badge](https://img.shields.io/static/v1?style=for-the-badge&label=GitHub&message=06&color=B48EAD&logo=github)](https://github.com/QCBSRworkshops/workshop06)

---

.center[
**Membres du CSBQ qui ont contribué à cet atelier**

en modifiant et améliorant son contenu dans le cadre du prix d’apprentissage et de développement
]

.pull-left[
.right[

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

[Laurie Maynard](https://github.com/LaurieMaynard)
 
[Pedro Henrique P. Braga](https://github.com/pedrohbraga)
 
[Katherine Hébert](https://github.com/katherinehebert)
 
[Alex Arkilanian](https://github.com/aarkilanian)
 
[Mathieu Vaillancourt](mathieu.vaillancourt.2@ulaval.ca)
 
[Esteban Góngora](https://github.com/estebangongora)

]
]

.pull-right[
.left[
**2019** - **2018** - **2017**

[Azenor Bideault](https://github.com/Azenor)
 
[Willian Vieira](https://github.com/willvieira)
 
[Pedro Henrique P. Braga](https://github.com/pedrohbraga)
 
[Marie Hélène Brice](https://github.com/mhBrice)
 
[Kevin Cazelles](https://github.com/KevCaz)

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

[Cédric Frenette Dussault]()
 
[Thomas Lamy]()
 
[Zofia Taranu](https://github.com/zetaranu)
 
[Vincent Fugère](https://github.com/VFugere)

]

]

.center[
__Si vous voulez contribuer aussi__, visitez [r.qcbs.ca/fr/contributing](https://r.qcbs.ca/fr/contributing/) et n'hésitez pas à [nous contacter](mailto:csbq.qcbs.r@gmail.com)!
]

---

# Objectifs d'apprentissage

1. Faire la distinction entre les modèles linéaires généralisés et les modèles généraux linéaires (incluant plusieurs de leurs équations!)

2. Identifier les situations où il est approprié d'utiliser des modèles linéaires généralisés.

3. Tester les hypothèses des modèles linéaires généralisés.

4. Implémenter et executer des modèles linéaires généralisés avec des données binaires, de proportion, et d'abondance.

5. Valider, interpreter and visualiser les résultats de modèles linéaires généralisés.

.small[
.center[
.alert[Important]
 
Cet atelier développera vos connaissances de modèles, y compris certaines équations importantes. *Ne vous sauvez pas!* Nous aborderons tout ceci progressivement, et vous sortirez de l'atelier avec plus de confiance!
]
]

---

# Packages requis

Pour suivre cet atelier, vous devriez avoir installé [R](https://cran.rstudio.com/) et [RStudio](https://rstudio.com/products/rstudio/download/#download).

.pull-left[
Cet atelier requiert aussi les packages R suivants:
1. [ggplot2](https://cran.r-project.org/package=ggplot2)
3. [MASS](https://cran.r-project.org/package=MASS)
4. [vcdExtra](https://cran.r-project.org/package=vcdExtra)
5. [bbmle](https://cran.r-project.org/package=bbmle)
6. [DescTools](https://cran.r-project.org/package=DescTools)
]

.pull-right[

Pour les installer à partir de CRAN, roulez:

```r
install.packages(c("ggplot2", 
                   'MASS', 
                   'vcdExtra',
                   'bbmle',
                   'DescTools')
                 )
```
]

Nous utiliserons aussi ces jeux de données:

- <a href="https://r.qcbs.ca/workshop06/pres-fr/data/mites.csv">mites.csv</a> 
- <a href="https://r.qcbs.ca/workshop06/pres-fr/data/faramea.csv">faramea.csv</a>

---
class: inverse, center, middle

# Préface

### Révision de modèles linéaires

---
# Répondre à des questions de recherche avec la modélisation

La plupart de nos recherches tentent d'expliquer des tendances dans nos observations à l'aide de variables prédictives.

Nous cherchons souvent une fonction `$f$` qui explique une variable réponse ( `$Y$` ) *en fonction* d'une ( `$X_1$` ) ou de plusieurs variables prédictives ( `$X_2$`, `$X_3$`, `$...$`, `$X_n$` ):
`$$Y = f(X_1)$$`
--
L'ensemble de variables prédictives que nous avons mesuré ne pourra jamais complètement expliquer notre variable `$Y$`. Il y a une **variation imprévisible** dans nos modèles, i.e. l'erreur `$\epsilon$`, qui fera toujours partie de notre fonction:

`$$Y = f(X_1, \epsilon)$$`
--
exclude:true
.small[
Dans l'[atelier 4](https://qcbsrworkshops.github.io/workshop04/pres-fr/workshop04-pres-fr.html#1), nous avons appris comment utiliser ces modèles linéaires généraux pour décrire la relation entre variables:

.center[
.pull-left[
.pull-left[
Test de `$t$`
]
.pull-right[
ANOVA
]
]
.pull-right[
.pull-left[
Régression linéaire
]
.pull-right[
ANCOVA
]
]
]
]
---
# Révision des modèles linéaires

Dans l'[atelier 4](https://qcbsrworkshops.github.io/workshop04/pres-fr/workshop04-pres-fr.html#1), nous avons appris comment utiliser les modèles linéaires généraux pour décrire la relation entre variables.

La formule générale de notre fonction linéaire `$Y = f(X_1)$` serait représentée comme:

`$$Y = \beta_0 + \beta_1X_i + \varepsilon_i$$`
où

`$Y_i$` est la valeur prédite de la variable réponse

`$\beta_0$` est le *coefficient inconnu* de l'**ordonnée à l'origine**

`$\beta_1$` est le *coefficient inconnu* de la **pente**

`$X_i$` est la valeur de la variable explicative

`$\varepsilon_i$` représente les résidus du modèle obtenus d'une distribution normale de moyenne 0 et de variance constante (qui est à estimer).

---

# Révision des modèles linéaires et leurs conditions

Nous avons aussi appris que les **modèles linéaires** produisent seulement des estimateurs non-biaisés (c'est-à-dire, sont seulement fiables) si ils suivent quelques conditions. Notamment:

1\. La population peut être décrite par une relation linéaire:

`$$Y = \beta_0 + \beta_1X_i + \varepsilon$$`
--

2\. Le terme d'erreur `$\varepsilon$` a la même variance quelque soit la valeur de la variable explicative (c'est-à-dire, l'homoscédasticité), et les termes d'erreur ne sont pas corrélés entre les observations (donc, il n'y a pas d'autocorrélation).

.pull-left[
`$$\mathbb{V}{\rm ar} (\epsilon_i | \mathbf{X} ) = \sigma^2_\epsilon,\ \forall i = 1,..,N$$`
]

.pull-right[
`$$\mathbb{C}{\rm ov} (\epsilon_i, \epsilon_j) = 0,\ i \neq j$$`
]

3\. Les résidus suivent une distribution normale:

`$$\boldsymbol{\varepsilon} | \mathbf{X} \sim \mathcal{N} \left( \mathbf{0}, \sigma^2_\epsilon \mathbf{I} \right)$$`

---
 
.small[
.center[**POURQUOI TOUTES CES ÉQUATIONS?**]
.center[*Courage! Ne perdez pas l'espoir! Elles sont nécessaires et nous expliquerons pourquoi!*]
]

---
## Qu'est-ce qu'on veut dire par:

.pull-left[
1\. La population peut être décrite par une relation linéaire:

`$$Y = \beta_0 + \beta_1X_i + \varepsilon$$`

2\.  Homoscédasticité et aucune autocorrélation

`$$\mathbb{V}{\rm ar} (\epsilon_i | \mathbf{X} ) = \sigma^2_\epsilon,\ \forall i = 1,..,N$$`
`$$\mathbb{C}{\rm ov} (\epsilon_i, \epsilon_j) = 0,\ i \neq j$$`

3\. Et, les résidus suivent une distribution normale:

`$$\boldsymbol{\varepsilon} | \mathbf{X} \sim \mathcal{N} \left( \mathbf{0}, \sigma^2_\epsilon \mathbf{I} \right)$$`

]

.pull-right[

Les estimations d'un modèle général linéaire telles que `$\widehat{Y} = \widehat{\beta}_0 + \widehat{\beta}_1 X$` assumemt que les données sont générées selon les conditions présentées.

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

???

Décrivez la figure aux participants, en montrant la ligne ajustée et les résidus, et en expliquant que la variance n'augmente pas avec les deux variables.

---
## Representer les conditions d'application

Simulons 250 unités d'observation qui satisfaient nos conditions d'application:  `$\epsilon_i \sim \mathcal{N}(0, 2^2), i = 1,...,250$`.

.pull-left[

```r
nSamples <- 250
ID <- factor(c(seq(1:nSamples)))

PredVar <- runif(nSamples, 
 min = 0, 
 max = 50)

simNormData <- data.frame(
 ID = ID,
 PredVar = PredVar,
 RespVar = (2*PredVar + 
 rnorm(nSamples,
 mean = 0,
 sd = 2)
 )
 )

# Nous avons appris à utiliser lm()

lm.simNormData <- lm(RespVar ~ PredVar, 
 data = simNormData)
```

]

.pull-right[

```r
layout(matrix(c(1,2,3,4),2,2)) 
plot(lm.simNormData)
```

]

???

Résidus: résidu =  y observé   –   y prédit par la modèle

Données aberrantes: observations avec des résidus larges, i.e. la valeur observée pour un point est très différente de la valeur prédite par le modèle linéaire.

Point de levier: un point de levier est défini comme une observation qui a une valeur de X qui est très éloignée de la moyenne de X.

Observations influentes: Une observation influente est défini par une observation qui change la pente de la relation. Donc, un point influent va avoir une influence élevée sur la prédiction du modèle. Une méthode pour trouver des observations influentes est de comparer un modèle avec et sans la dite observation.

1. Ces graphiques permettent de vérifier les conditions d'application de la linéarité et l'homoscédasticité.

2. Le graphique QQ permet la comparaison des résidus avec une distribution normal.

3. Scale-location plot (square rooted standardized residual vs. predicted value) est utile pour vérifier l'homoscédasticité;

4. La distance de Cook est une mesure d'influence des observations sur le coefficient de la régression linéaire.

---
# Un exemple de modèle linéaire général

Utilisons nos connaissances des modèles linéaires généraux pour explorer la relation entre les variables dans le *jeu de données de mites Orbatid*.

.pull-left3[
Commençons par charger les données:

.small[

```r
mites <- read.csv('data/mites.csv', 
 stringsAsFactors = TRUE)
```

]
]

.pull-right3[
Et ensuite, les explorer:

```r
str(mites)
# essayez aussi head(mites)
```
]

--

.small[

```
# 'data.frame':	70 obs. of  9 variables:
#  $ Galumna   : int  8 3 1 1 2 1 1 1 2 5 ...
#  $ pa        : int  1 1 1 1 1 1 1 1 1 1 ...
#  $ totalabund: int  140 268 186 286 199 209 162 126 123 166 ...
#  $ prop      : num  0.05714 0.01119 0.00538 0.0035 0.01005 ...
#  $ SubsDens  : num  39.2 55 46.1 48.2 23.6 ...
#  $ WatrCont  : num  350 435 372 360 204 ...
#  $ Substrate : Factor w/ 7 levels "Barepeat","Interface",..: 4 3 2 4 4 4 4 2 3 4 ...
#  $ Shrub     : Factor w/ 3 levels "Few","Many","None": 1 1 1 1 1 1 1 2 2 2 ...
#  $ Topo      : Factor w/ 2 levels "Blanket","Hummock": 2 2 2 2 2 2 2 1 1 2 ...
```

]

.small[

Le jeu de données `mites` contient **70 échantillons de mousses et mites**, **5 variables environmentales**, l'abondance de la mite *Galumna* sp., et l'abondance totale des mites.

]

???

Il s'agit d'une partie du jeu de données classique des
[mites Oribatid (Acari,Oribatei)](http://adn.biol.umontreal.ca/~numericalecology/data/oribates.html) des sphaignes (*Sphagnum* sp.) du Lac Geai, [Station de Biologie de l'Université de Montréal](https://goo.gl/maps/PxN1Q7KUPnUt92Eu5).
]

---
# Quelles questions pouvons-nous poser?

**Question**: *Est-ce que l'environnement permet de prédire l'abondance, l'occurrence, ou la proportion de Galumna sp.?*

.pull-left[
.small[
Variables réponses:
1. Occurrence: `pa` 
2. Abondance: `Galumna`
3. Proportion ou fréquence relative: `prop`
]
]
.pull-right[
.small[
Variables prédictives:
1. Densité du substrat: `SubsDens`
2. Contenu en eau (du sol): `WatrCont`
3. Substrat: `Substrate`
4. Arbustes: `Shrub`
5. Topographie: `Topo`
]
]

**Pour répondre à notre question**: *On peut développer des modèles!*

.pull-left[
`$\text{Abondance} = f(\text{Contenu en eau}, \epsilon)$`

`$\text{Proportion} = f(\text{Contenu en eau}, \epsilon)$`

`$\text{Occurrence} = f(\text{Substrat}, \epsilon)$`

`$\text{Abondance} = f(\text{Topographie}, \epsilon)$`
]
.pull-right[
`$\text{Occurrence} = f(\text{Arbustes}, \epsilon)$`

`$\text{Fréquence relative} = f(\text{Topographie}, \epsilon)$`

`$\text{Occurrence} = f(\text{Densité du substrat}, \epsilon)$`

`$\text{Abondance} = f(\text{Substrat}, \epsilon)$`
]

---
exclude: true

# Par où commencer?

.small[
Pouvons-nous voir une relation entre *Galumna* et une ou plusieurs des cinq variables environnementales?
]

.pull-left2[

]
.pull-right2[

```r
plot(mites)
```

Peut-être...

`Galumna` `$vs$` `WatrCont`?!
]

---
# Explorer les relations

Est-ce que la communauté de _Galumna_ (abondance, occurence, et fréquence relative) varie en fonction du contenu en eau du sol?

.small[

```r
par(mfrow = c(1, 3), cex = 1.4)
plot(Galumna ~ WatrCont, data = mites, 
     xlab = 'Contenu en eau', ylab = 'Abondance')
boxplot(WatrCont ~ pa, data = mites, 
        xlab = 'Presence-absence', ylab = 'Contenu en eau', col = 'red')
plot(prop ~ WatrCont, data = mites, 
     xlab = 'Contenu en eau', ylab = 'Frquence relative')
```

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

---

# Utiliser des modèles linéaires

Utilisons des modèles linéaires (généraux) pour tester si `Galumna`, `pa`, et/ou `prop`, varient en fonction du contenu d'eau (`WatrCont`) en utilisant la fonction `lm()`:

--
.small[
.pull-left3[

```r
lm.abund <- lm(Galumna ~ WatrCont, data = mites)
## summary(lm.abund)

lm.pa <- lm(pa ~ WatrCont, data = mites)
## summary(lm.pa)

lm.prop <- lm(prop ~ WatrCont, data = mites)
## summary(lm.prop)
```
]

.pull-right3[
`$\text{Galumna} = f(\text{Contenu en eau}, \epsilon)$`

`$\text{Occurrence} = f(\text{Contenu en eau}, \epsilon)$`

`$\text{Proportion} = f(\text{Contenu en eau}, \epsilon)$`
]
]

.small[
.pull-left3[

```r
# Extraction du Pr(>|t|)

summary(lm.abund)$coefficients[, 4]
# (Intercept) WatrCont 
# 3.981563e-08 1.206117e-05
summary(lm.pa)$coefficients[, 4]
# (Intercept) WatrCont 
# 6.030252e-12 4.676755e-08
summary(lm.prop)$coefficients[, 4]
# (Intercept) WatrCont 
# 4.977432e-08 1.665437e-05
```
]
]
.center[
.pull-right3[
 
Tous les modèles sont significatifs!

.alert[Mais...] **avons-nous respecté les conditions du modèle?**
]
]

---

# Utiliser des modèles linéaires: conditions

.alert[Attends une minute... on a un problème!]:

.small[
.pull-left[

```r
plot(Galumna ~ WatrCont, data = mites)
abline(lm.abund)
```

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

```r
par(mfrow = c(2, 2), cex = 1.4)
plot(lm.abund)
```

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

---

# Utiliser des modèles linéaires: conditions

Encore pire pour les autres modèles (`lm(prop ~ WatrCont)`) :

.pull-left[

```r
plot(prop ~ WatrCont, data = mites)
abline(lm.prop)
```

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

.pull-right[

```r
par(mfrow = c(2, 2), cex = 1.4)
plot(lm.prop)
```

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

---

# Utiliser des modèles linéaires: conditions

Encore pire pour les autres modèles (`lm(pa ~ WatrCont)`) :

.pull-left[

```r
plot(pa ~ WatrCont, data = mites)
abline(lm.pa)
```

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

.pull-right[

```r
par(mfrow = c(2, 2), cex = 1.4)
plot(lm.pa)
```

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

---
# Révision: Conditions d'application du modèle linéaire

.small[
.pull-left[

Rappelez-vous de notre modèle linéaire de base?

`$$Y_i = \beta_0 + \beta_1X_i + \varepsilon$$`
On peut aussi l'écrire comme:

`$$Y_i \sim N(\mu = \beta_0 + \beta_1 X_i +\varepsilon, \sigma^2)$$`

avec `$N(\cdot)$`, c'est-à-dire que `$Y_i$` est échantillonné d'une **distribution normale**avec les paramètres `$\mu$` (moyenne; qui dépend de `$X_i$`) et `$\sigma$` (variance; qui a la même valeur pour tous les `$Y_i$`).
]

.pull-right[
3 valeurs pour `$\mu$`, `$\sigma = 5$`

`$\mu = 25$`, 3 valeurs pour `$\sigma$`
<img src="workshop06-pres-fr_files/figure-html/unnamed-chunk-21-1.png" width="216" style="display: block; margin: auto;" />
]
]
???
Sur la droite, nous montrons des figures où mu varie ou oméga varie. La figure du haut marcherait avec un modèle linéaire puisque mu varie, mais pas oméga. La figure du bas montre qu'on ne peut utiliser un modèle linéaire puisque oméga varie. Il faudrait donc explorer dse distribtuions différentes. 
---
# Prédiction du modèle

Rappel: nous voulons estimer les *coefficients inconnus* `$\beta_0$` et `$\beta_1$`, pour tracer une ligne droite qui prédit chaque valeur de `$Y$` en fonction de `$X$`!

`$$Y_i \sim N(\mu = \beta_0 + \beta_1 X_i +\varepsilon, \sigma^2)$$`

On peut obtenir les paramètres `$\mu$` et `$\sigma^2$` pour une distribution normale correspond à notre équation!

Pour obtenir les coefficient de nos modèles, on peut utiliser la fonction `coef()`:

.pull-left[

```r
coef(lm.abund)
#  (Intercept)     WatrCont 
#  3.439348672 -0.006044788
```
]

.pull-right[

```r
summary(lm.abund)$sigma
# [1] 1.513531
```
]

Quels sont les paramètres de la distribution normale utilisée pour modéliser `$Y$` quand le contenu en eau est `$300$`?

.center[
pour `$X_{i} = 300$`, `$\mu = 3.44 + (-0.006 \times 300) = 1.63$`   et    `$\sigma^2 = 1.51$`
]

c'est-à-dire, des valeurs `$Y$` échantillonés aléatoirement quand le contenu en eau est de 300 devraient avoir une moyenne de `$1.63$` et une variance de `$1.51$`.

---
# Prédiction du modèle

.pull-left[
`$$y_i \sim N(\mu = \beta_0 + \beta_1 X_i +\varepsilon, \sigma^2)$$`
`$$\mu = 3.44 + (-0.006 \times 400) +\varepsilon= 1.02$$`
`$$\sigma^2 = 1.51$$`
]
.pull-right[
Lorsque `$x_i = 300$`, `$Y_i$` suit une distribution normale avec `$\mu = 1.63$` et `$\sigma^2 = 1.51$`.

Lorsque `$x_i = 400$`, `$Y_i$` suit une distribution normale avec `$\mu = 1.02$` et `$\sigma^2 = 1.51$`, etc.
]

.pull-left[

.small[Notre modèle prévoit donc que les observations se situent dans les zones grisées suivantes:]

.center[
<img src="workshop06-pres-fr_files/figure-html/unnamed-chunk-22-1.png" width="324" style="display: block; margin: auto;" />
]
]

.pull-right[

.center[**Mais, le font-elles ?**]

]

---
# Prédiction du modèle

Lorsque `$x_i = 400$`, `$Y_i$` suit une distribution normale avec `$\mu = 1.02$` et `$\sigma^2 = 1.51$`, etc.
]

.pull-left[

.small[Notre modèle prévoit donc que les observations se situent dans les zones grisées suivantes:]

.center[
<img src="workshop06-pres-fr_files/figure-html/unnamed-chunk-23-1.png" width="324" style="display: block; margin: auto;" />

]
]

.pull-right[

.center[**Mais, le font-elles ?**]

.center[*Hmmmmmmmm...*]

**Lesquelles des propositions suivantes sont correctes?**

1. Les résidus semblent violer les premises des modèles linéaires ;
2. La `$\sigma^2$` (variance) change selon les valeurs de `$x$` ;
3. Les observations ne tombent pas sur les valeurs attendues.
]

???

Vous pouvez effectuer un sondage à choix multiples auprès des participants.

Les options 1, 2 et 3 posent problème dans ce modèle, c'est-à-dire qu'elles sont toutes correctes.

---
## Comment règler ces problèmes? Transformer?

.small[
Souvent, les données ne respectent pas les conditions nous venons de voir, montrant des traces de **non-normalité** et/ou d'**hétéroscédasticité**.

On se fait souvent dire qu'il faut **transformer** nos données en utilisant des transformations logarithmiques, racine carrée, et cosinus.
]

.pull-left[
Mais, les transformations ne règlent pas toujours le problème, et viennent avec quelques inconvénients:

1. Elles changent la variable réponse, ce qui peut compliquer l'interprétation;  
2. Elles ne peuvent pas toujours améliorer la linéarité et l'homogénité de la variance en même temps;  
3. Les limites de l'espace d'échantillonage changent.

]

.pull-right[
Par exemple, notre modèle linéaire:

`$$Y_i = \beta_0 + \beta_1X_i + \varepsilon$$`
devient ceci après une transformation logarithmique:

$$E(\log{Y_i}) = \beta_0 + \beta_1X_i $$

*C'est moins intuitif d'interpréter que l'abondance de Galumna prend la forme de `$\log(1.63).$` pour chaque augmentation de `$300$` unités en contenu en eau!*
]

---
# En plus, la distribution normale n'est pas la seule option!

Les statisticiens ont développé [de nombreuses lois de probabilité (distributions)](https://www.causascientia.org/math_stat/Dists/Compendium.pdf) pour décrire divers types de données.

Une distribution donne la probabilité d'observer chaque issue possible d'une expérience ou échantillonage (e.g. `$abund = 8$` individus de *Galumna*).

Les distributions peuvent être **discrètes** (que des nombres entiers), **continues** (incluent aussi des fractions) et bornées par une certaine gamme de valeurs (e.g., `$0$` à `$1$`, `$0$` à `$+\infty$`, `$-\infty$` à `$+\infty$`).

Toutes les distributions ont des **paramètres** qui déterminent leur forme (e.g. `$\mu$` et `$\sigma^2$` pour la distribution normale).

---
# Données biologiques & distributions: Poisson

L'abondance de *Galumna* suit une distribution discrète (que des nombres entiers).

Pour modéliser les données d'abondance, la [distribution de Poisson](https://fr.wikipedia.org/wiki/Loi_de_Poisson) est souvent utilisée:

- **Poisson** est une distribution discrète avec un seul paramètre, `$\lambda$` (lambda), qui détermine la moyenne et la variance de la distribution:

---
# Données biologiques & distributions: Poisson

*Galumna* semble suivre une loi de Poisson avec une faible valeur de `$\lambda$` :

.pull-left[

```r
  hist(mites$Galumna)
```

<img src="workshop06-pres-fr_files/figure-html/unnamed-chunk-25-1.png" width="432" style="display: block; margin: auto;" />
]
.pull-right[
Est-ce qu'on modèle avec une distribution de Poisson fonctionnerait mieux pour le modèle `$\text{Galumna} = f(\text{Water content})$` ?
<img src="workshop06-pres-fr_files/figure-html/unnamed-chunk-26-1.png" width="432" style="display: block; margin: auto;" />
]

---
# Données biologiques & distributions: Poisson

Les données de présence-absence prennent une autre forme :

- Seulement des `0`s et des `1`s
- La loi de Poisson n'est pas appropriée pour cette variable

```r
  hist(mites$pa)
```

---
# Données biologiques & distributions: Bernouilli

**Distribution de [Bernoulli](https://fr.wikipedia.org/wiki/Loi_de_Bernoulli)** :

- N'inclut que deux issues possibles dans son ensemble: succès (`1`) ou échec (`0`)
- N'a qu'un paramètre, `$p$`, la probabilité de succès

Nous pouvons utiliser la loi de Bernoulli pour calculer la probabilité d'obtenir l'issue Galumna "present" (`1`) vs. "absent" (`0`).

---
# Données biologiques & distributions: Binomiale

**[Distribution binomiale](https://fr.wikipedia.org/wiki/Loi_binomiale)** : Lorsqu'il y a plusieurs épreuves (chacune avec un succès/échec), la loi de Bernoulli devient la loi binomiale

- Inclut le paramètre additionel `$n$`, le nombre d'épreuves
- Prédit la probabilité d'observer une certaine proportion de succès, `$p$`, sur le nombre total d'épreuves, `$n$`

---
# Données biologiques & distributions: Binomiale

**Distribution binomiale** : utilisée pour modéliser des données lorsque le nombre de succès est donné par un nombre entier, et lorsque le nombre d'épreuves, `$n$`, est connu.

**Différence principale avec la loi de Poisson** : L'étendue de la loi binomiale a une limite supérieure, `$n$`. Par conséquent, elle est asymétrique et décalée à gauche lorsque `$p$` est faible, mais décalée à droite lorsque `$p$` est élevé.

---
class: middle, center

*Retournons à notre problème... que faire si nous avons d'autres distributions?*

.center[On peut utiliser un **modèle linéaire généralisé**!]

Dans`R`, on va donc passer de la fonction pour les modèles linéaires (généraux) `lm()` et commencerons à utiliser la fonction `glm()`!

---
# Modèles linéaires généralisés

&nbsp; Un **modèle linéaire généralisé** est composé d'un **prédicteur linéaire**:

`$$\underbrace{g(\mu_i)}_{Fonction~~lien}  = \underbrace{\beta_0 + \beta_1X_1~+~...~+~\beta_pX_p}_{Composante~linéaire}~+~\underbrace{\epsilon}_{Error}$$`
&nbsp; et,

&nbsp; 1. une fonction **lien** `$g(\mu_i)$` qui **transforme la valeur prédite** de `$Y$` et décrit comment la moyenne `$E(Y_i) = \mu_i$` dépend du prédicteur linéaire

`$$g(\mu_i) = ηi$$`

&nbsp; 2\. une fonction de **variance** qui décrit comment la variance `$\text{var}(Y_i)$` dépend de la moyenne

`$$\text{var}(Y_i) = \phi V(\mu)$$`
&nbsp; où le **paramètre de dispersion** `$phi$` est une constante.

???
La fonction de lien met en relation la valeur attendue de la variable réponse au prédicteur linéaire du modèle, pour que la relation entre le prédicteur et la réponse soit modélisée avec une régression linéaire.

Comment cela se différencie de transformer une variable?

La transformation la plus commune est une transformation en log puisque la valeur de y n'a pas l'Apparence de suivre une distribution normale.

Lorsque nous transformons les données d'un modèle linéaire général, nous n'affirmons plus que y est distribuer normalement autour de la moyenne pour une valeur de x donnée. Nous affirmons en fait que la nouvelle variable (log(y) ou ln(y)) est ditribuée normalement.

Dans le cas d'un modèle linéaire généralisé de Poisson, cependant, la fonction de lien ne change pas la distribution des observations actuelles pour tenter de les rendre normale. Plutôt, la fonction lien redéfinie la relation de la variable X directement avec la moyenne de y sous une distribution de Poisson. Les observations individuelles vont donc varier autour de la valeur attendue.

La moyenne du log n'est pas la le log de la moyenne! Vous ne pouvez prendre l'exposant de la moyenne de log(y) ou ln(y) et avoir la moyenne de y. D'un autre côté, la fonction lien vous permet de le faire!

---
# Le modèle linéaire général en GLM

Alors, pour notre **modèle linéaire général** avec `$\epsilon ∼ N(0, σ^2)$` et `$Y_i \sim{} N(\mu_i, \sigma^2)$`, on a le même prédicteur linéaire:

`$$\underbrace{g(\mu_i)}_{Fonction~~lien}  = \underbrace{\beta_0 + \beta_1X_1~+~...~+~\beta_pX_p}_{Composante~linéaire}~+~\underbrace{\epsilon}_{Error}$$`
la fonction lien: 
`$$g(\mu_i) = \mu_i$$`

et, la fonction de variance:

`$$V(\mu_i) = 1$$`

.center[*Wow! C'est équivalent au modèle linéaire normal que nous avions vu au début de l'atelier!*]

---

# `glm()`

Dans `R`, on peut estimer un modèle linéaire généralisé avec la fonction `glm()`, qui ressemble beaucoup à la fonction `lm()`: 
.pull-left2[

```r
glm(formula, 
    family = gaussian(link = "identity"), 
    data,
    ...)
```
]
.pull-right2[

```r
lm(formula, 
    data,
    ...)
#
```

]

avec l'argument `family` (voir `?family`) qui prend le nom de la **fonction lien** (à l'intérieur de `link()`) et la **fonction de variance**.

Cette approche s'applique à d'autres distributions!

???
La fonction d'identité pour une modèle linéaire général est 'identité'(voir la diapo suivante)

| Distribution de `$Y$` | Nom du fonction lien | Fonction lien | Modèle      | `R` |
|---------------------|--------------------|---------------|------------|-----|
| Normale              | Identité           | `$g(\mu) = \mu$`  | `$\mu = \mathbf{X} \boldsymbol{\beta}$` | `gaussian(link="identity")` | 
| Binomiale            | Logit              | `$g(\mu) = \log\left(\dfrac{\mu}{1-\mu}\right)$`  | `$\log\left(\dfrac{\mu}{1-\mu}\right) = \mathbf{X} \boldsymbol{\beta}$` | `binomial(link="logit")`| 
| Poisson              | Log           | `$g(\mu) = \log(\mu)$` | `$\log(\mu) = \mathbf{X} \boldsymbol{\beta}$` | `poisson(link="log")`|
| Exponentielle              | Inverse négative          | `$g(\mu) = -\mu^{-1}$`  | `$-\mu^{-1} = \mathbf{X} \boldsymbol{\beta}$` | `Gamma(link="inverse")` |

???

Comment cela se différencie de transformer une variable?

La transformation la plus commune est une transformation en log puisque la valeur de y n'a pas l'Apparence de suivre une distribution normale.

La moeynne du log n'est pas la le log de la moyenne! Vous ne pouvez prendre l'exposant de la moyenne de log(y) ou ln(y) et avoir la moyenne de y. D'un autre côté, la fonction lien vous permet de le faire!
---
class: inverse, center, middle

# GLM with avec des données binaire
## Modèles linéaires généralisés

---
# Variables binaires

Une variable réponse commune dans les jeux de données en écologie est la variable binaire.

On observe un phénomène `$Y$` ou son "absence":

- Présence/Absence d'une espèce
- Présence/Absence d'une maladie
- Succès/Échec d'observer un comportement
- Survie/Mort d'un organisme

Souvent, on veut savoir si l'occurrence d'espèces varient en fonction de l'environnement

`$$\text{Occurrences} = f(\text{Environment})$$`

---
# Variables binaires

Dans un modèle linéaire, les valeurs prédites peuvent se trouver hors de l'intervalle `[0, 1]` avec `lm()`:

```r
lm(Pres ~ ExpVar)
```

.pull-left[
<img src="workshop06-pres-fr_files/figure-html/unnamed-chunk-34-1.png" width="540" style="display: block; margin: auto;" />
]

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

---
# GLM avec données binomiales: lien logit

Un **modèle linéaire généralisé** est composé d'un **prédicteur linéaire**:

`$$\underbrace{g(\mu_i)}_{Fonction~~lien}  = \underbrace{\beta_0 + \beta_1X_1~+~...~+~\beta_pX_p}_{Composante~linéaire}~+~\underbrace{\epsilon}_{Error}$$`

Pour des données binaires ou binomiales, la fonction lien est nommée le **logit**:

`$$g(p) = log\frac{p}{1-p}$$`

où `$\mu$` est la valeur attendue (probabilité que `$Y = 1$`).

La fonction *logit* est le logarithme de la probabilité `$(\frac{p}{1-p})$`.

Cette probabilité place nos valeurs attendues entre `0` et `+Inf`.

La transformation log place nos valeurs attendues entre `-Inf` et `+Inf`.

*Les valeurs attendues sont maintenant liées de façon **linéaire** au prédicteur linéaire.*

???

On peut comprendre que si l'évènement A a une probabilité p de se produire, la chance que l'évènement A se produise est le ratio de la probabilité que A se produise sur la probabilité que A ne se produise pas: p/(1−p).

Par exemple, si j'ai une probabilité de 0.6 pour échouer à mon cours, la chance que j'échouerai mon cours est 0.6/(1 − 0.6) = 1.5. Donc, la probabilité d'observer un échec est 1.5 fois plus grand que la probabilité de ne pas en observer (donc, 1.5 × 0.4 = 0.6).

---
# GLM avec données binomiales: lien logit

Un rappel de la fonction `glm()`:

```r
glm(formula, 
    family = ???, 
    data,
    ...)
```

| Distribution de `$Y$` | Nom du fonction lien | Fonction lien | Modèle      | `R` |
|---------------------|--------------------|---------------|------------|-----|
| Normale              | Identité           | `$g(\mu) = \mu$`  | `$\mu = \mathbf{X} \boldsymbol{\beta}$` | `gaussian(link="identity")` | 
| Binomiale            | Logit              | `$g(\mu) = \log\left(\dfrac{\mu}{1-\mu}\right)$`  | `$\log\left(\dfrac{\mu}{1-\mu}\right) = \mathbf{X} \boldsymbol{\beta}$` | `binomial(link="logit")`| 
| Poisson              | Log           | `$g(\mu) = \log(\mu)$` | `$-\mu^{-1} = \mathbf{X} \boldsymbol{\beta}$` | `poisson(link="log")`|
| Exponentielle              | Inverse négative          | `$g(\mu) = -\mu^{-1}$`  | `$\log(\mu) = \mathbf{X} \boldsymbol{\beta}$` | `Gamma(link="inverse")` |

```r
glm(formula, 
    family = binomial(link = "logit"), # aussi nommé "logistique"
    data,
    ...)
```

---
# Exercice 1 - Notre premier modèle linéaire, ensemble.

Spécifiez un modèle de régression logistique avec le jeu de données `mites`.

```r
#setwd('...')
mites <- read.csv("mites.csv", header = TRUE)
str(mites)
```

???

Un rappel du jeu de données:

70 sites (rangés) x 35 morphospecies (colones). Les unités d'échantillonage sont des coeurs de substrat de 5 cm en diamètre et 10 cmen profondeur.

SubsDens: densité du substratum en g.L-1 de matière sèche décompressée (variable quantitative)

WatrCont: contenu en eau du substratum, en pourcentage de volume (variable quantitative)

Borcard, D., P. Legendre and P. Drapeau. 1992. Partialling out the spatial component of ecological variation. Ecology 73: 1045-1055.

Borcard, D. and P. Legendre. 1994. Environmental control and spatial structure in ecological communities: an example using Oribatid mites (Acari, Oribatei). Environmental and Ecological Statistics 1: 37-53.

---
# Exercice 1

Spécifions un modèle de la présence et de l'absence de *Galumna sp.* en fonction du contenu en eau du sol et de la topographie.

```r
logit.reg <- glm(pa ~ WatrCont + Topo, data=mites,
 family = binomial(link = "logit"))
```

```r
summary(logit.reg)
```

???

Mettez de l'emphase sur l'argument `family`,  et référez à ?family.

---
# Exercice 1

.small[

```r
summary(logit.reg)
# 
# Call:
# glm(formula = pa ~ WatrCont + Topo, family = binomial(link = "logit"), 
#     data = mites)
# 
# Deviance Residuals: 
#     Min       1Q   Median       3Q      Max  
# -2.0387  -0.5589  -0.1594   0.4112   2.0252  
# 
# Coefficients:
#              Estimate Std. Error z value Pr(>|z|)    
# (Intercept)  4.464402   1.670622   2.672 0.007533 ** 
# WatrCont    -0.015813   0.004535  -3.487 0.000489 ***
# TopoHummock  2.090757   0.735348   2.843 0.004466 ** 
# ---
# Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# 
# (Dispersion parameter for binomial family taken to be 1)
# 
#     Null deviance: 91.246  on 69  degrees of freedom
# Residual deviance: 48.762  on 67  degrees of freedom
# AIC: 54.762
# 
# Number of Fisher Scoring iterations: 6
```
]

???

En expliquant cette sortie, montrez ce que sont les blocs et rappelez-vous que sa structure ressemble à celle de `summary.lm()`, mais qu'il y a quelques différences particulières (par exemple le paramètre de dispersion).

La déviation, les coefficients et la dispersion seront abordés dans les diapositives suivantes.

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

.small[
En utilisant le jeu de données `bacteria`, spécifiez un modèle de la présence de *H. influenzae* en fonction du traitement et de la semaine de test.

Commencez avec un modèle complet et trouvez le modèle le plus parcimonieux.

```r
library(MASS)
data(bacteria)
str(bacteria)
# 'data.frame':	220 obs. of  6 variables:
#  $ y   : Factor w/ 2 levels "n","y": 2 2 2 2 2 2 1 2 2 2 ...
#  $ ap  : Factor w/ 2 levels "a","p": 2 2 2 2 1 1 1 1 1 1 ...
#  $ hilo: Factor w/ 2 levels "hi","lo": 1 1 1 1 1 1 1 1 2 2 ...
#  $ week: int  0 2 4 11 0 2 6 11 0 2 ...
#  $ ID  : Factor w/ 50 levels "X01","X02","X03",..: 1 1 1 1 2 2 2 2 3 3 ...
#  $ trt : Factor w/ 3 levels "placebo","drug",..: 1 1 1 1 3 3 3 3 2 2 ...
```

*Tests de la présence de la bactérie H. influenzae dans les enfants avec une otite moyenne dans le Territoire du Nord de l'Australie*

]

???

Ceci est un exercise de groupe. Formez des groupes de 3-4 personnes dans des salles de réunion et donnez-leur 5 minutes pour écrire le code en collaboration.

Plus d'infos sur le jeu de données:

Dr A. Leach a testé les effets d'un un médicament sur 50 enfants ayant des antécédents d'otite moyenne dans le Territoire du Nord de l'Australie. Les enfants ont été choisis aléatoirement pour recevoir le médicament ou un placebo.

La présence de H. influenzae a été contrôlée aux semaines 0, 2, 4, 6 et 11 : 30 de ces contrôles étaient manquants et ne sont pas inclus dans ce jeu de données.

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

```r
model.bact1 <- glm(y ~ trt * week, data = bacteria, family = binomial)
```

```r
model.bact2 <- glm(y ~ trt + week, data = bacteria, family = binomial)
```

```r
model.bact3 <- glm(y ~ week, data = bacteria, family = binomial)
```

```r
anova(model.bact1, model.bact2, model.bact3, test = "LRT")
# Analysis of Deviance Table
# 
# Model 1: y ~ trt * week
# Model 2: y ~ trt + week
# Model 3: y ~ week
#   Resid. Df Resid. Dev Df Deviance Pr(>Chi)  
# 1       214     203.12                       
# 2       216     203.81 -2  -0.6854  0.70984  
# 3       218     210.91 -2  -7.1026  0.02869 *
# ---
# Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
```

???

LRT fait un test de ratio de likelihood (LR) entre modèles. On l'obtient dans l'exemple ci-dessous.

# Improve this in the future

vals <- (sum(residuals(base)^2) - sum(residuals(full)^2))/sum(residuals(full)^2) * full$df.residual

pchisq(vals, df.diff, lower.tail = FALSE)

---
# Interpréter les coefficients de `glm()`

Regardez à nouveau les coefficients du modèle `logit.reg` :

```r
logit.reg <- glm(pa ~ WatrCont + Topo, data=mites,
 family = binomial(link = "logit"))
```

```r
summary(logit.reg)$coefficients
#                Estimate  Std. Error   z value     Pr(>|z|)
# (Intercept)  4.46440199 1.670622482  2.672299 0.0075333598
# WatrCont    -0.01581255 0.004535069 -3.486728 0.0004889684
# TopoHummock  2.09075654 0.735348234  2.843220 0.0044660283
```

La sortie indique que le contenu d'eau et la topographie sont associés significativement avec la variable présence-absence.

.center[.comment[Mais comment interprète-on les coefficients de la pente?]]

---
# Interpréter les coefficients de `glm()`

L'interprétation directe des coefficients dans le modèle logit est délicate en raison de la fonction de lien. Si le lien est *identité*, l'interprétation est beaucoup plus facile.

Supposons que nous avons un résultat binaire `$y$` et deux covariables `$x_1$` et `$x_2$` (et une constante). La probabilité d'un résultat positif ( `$y = 1$` ) est donnée par:

`$$Pr(y_i = 1) = p = g^{-1}(\beta_0 + x_{1i}\beta_1 + x_{2i}\beta_2)$$`

où `$g^{-1}()$` est la **fonction de lien inverse**.

Maintenant, concentrons-nous sur l'interprétation du coefficient `$\beta_1$`.

#### Lien d'identité

Pour le lien d'identité, l'interprétation est simple. Pour l'augmentation d'une unité en `$x_1$`, `$\beta_1$` entraîne une différence constante dans le résultat.

`$\Delta{y_i} = [\beta_0 + (\color{red}{x_{1i} + 1})\beta_1 + x_{2i}\beta_2] - (\beta_0 + x_{1i}\beta_1 + x_{2i}\beta_2)$`

`$\Delta{y_i} = \beta_1$`

---
# Interpréter les coefficients de `glm()`

#### Lien logit

Si un modèle logistique linéaire utilise les deux covariables `$x_1$` et `$x_2$`, nous avons le modèle :

`$$log({\dfrac{p}{1-p}})=\beta_0 + x_{1i}\beta_1 + x_{2i}\beta_2$$`

pour une probabilité logarithmique de réponse positive. Nous pouvons également écrire ce modèle comme suit :

`$$\dfrac{p}{1-p}=exp(\beta_0 + x_{1i}\beta_1 + x_{2i}\beta_2)$$`

`$$Pr(yi) = \dfrac{exp(\beta_0 + x_{1i}\beta_1 + x_{2i}\beta_2)}{1 + exp{(\beta_0 + x_{1i}\beta_1 + x_{2i}\beta_2)}}$$`

.center[Super ! Lâche pas !]

---
# Interpréter les coefficients de `glm()`

#### Lien logit

Comme le lien inverse est non linéaire, l'interprétation du coefficient est difficile.

Voyons ce qui arrive aux coefficients pour un changement d'une unité à `$x_1$` :

`$$\Delta{y_i} = \dfrac{\exp[\beta_0 + (\color{red}{x_{1i} + 1})\beta_1 + x_{2i}\beta_2]}{1 + \exp{[\beta_0 + (\color{red}{x_{1i} + 1})\beta_1 + x_{2i}\beta_2]}} - \dfrac{\exp[\beta_0 + x_{1i}\beta_1 + x_{2i}\beta_2]}{1 + \exp{[\beta_0 + x_{1i}\beta_1 + x_{2i}\beta_2]}}$$`
$$\Delta{y_i} = exp(\beta_1) $$

Lorsque `$x$` augmente d'une unité, la probabilité augmente d'un facteur de `$\exp(\beta_1)$`.

---
# Interpréter les coefficients de `glm()`

###### Revenons à notre modèle `logit.reg`!

.small[

```r
logit.reg
# 
# Call:  glm(formula = pa ~ WatrCont + Topo, family = binomial(link = "logit"), 
#     data = mites)
# 
# Coefficients:
# (Intercept)     WatrCont  TopoHummock  
#     4.46440     -0.01581      2.09076  
# 
# Degrees of Freedom: 69 Total (i.e. Null);  67 Residual
# Null Deviance:	    91.25 
# Residual Deviance: 48.76 	AIC: 54.76
```
]

Pour une augmentation (ou une diminution) d'une unité de contenu en eau, on peut obtenir la probabilité de la présence de mites.

```r
exp(logit.reg$coefficient[2])
#  WatrCont 
# 0.9843118
```

---
# Interpréter les coefficients de `glm()`

.pull-left[

Quand _odds_ `$= 1$`, la probabilité d'observer l'événement `$Y$` est égale à la probabilité de *ne pas* l'observer (*i.e.* `$p = 0.5$`, donc `$0.5/(1-0.5) = 1$`).

Lorsque _odds_ `$< 1$`, il faut prendre l'inverse (_i._ `$1$` divisé par les _odds_).

L'interprétation est alors de savoir quelle est la **moindre** probabilité d'observer l'événement d'intérêt.

Pour la teneur en eau, le _odds_ est de `$0.984$`. L'inverse est :

`$$\dfrac{1}{0.984} = 1.0159$$`
_i.e._, il y a une augmentation d'une unité de la teneur en eau diminue la probabilité d'observer _Galumna sp_. de `$1.0159$`.

]

.pull-right[
 Nous pouvons obtenir une variation en pourcentage comme ci-dessous :
`$$(1.0159 - 1) * 100 = 1.59%$$`
Il y a une diminution de `$1.59%$` de la probabilité d'observer _Galumna sp._ avec une augmentation d'une unité de la teneur en eau.

]

---
# Pouvoir prédictif et ajustement du modèle

Rappelez-vous que la valeur `$R^2$` est une mesure du pouvoir explicatif du modèle.

Dans les modèles linéaires simples, `$R^2$` peut être obtenu comme le carré du coefficient de l'équation de corrélation de l'échantillon ($r$) entre les résultats observés et les valeurs prédictives.
--
 
 
Dans certains des modèles linéaires généralisés, nous adoptons d'autres approches.

Nous pouvons obtenir un [ `$\text{pseudo-R}^2$` ](https://stats.idre.ucla.edu/other/mult-pkg/faq/general/faq-what-are-pseudo-r-squareds/), un concept analogue au `$R^2$` pour les modèles estimés par maximisation de la vraisemblance.

Nous pouvons calculer le `$\text{pseudo-R}^2$` de McFadden (1973) :

`$$\text{pseudo-R}^2 = \frac{\text{déviance nulle - déviance résiduelle}}{\text{déviance nulle}}$$`

???
`$\text{pseudo-R}^2$` est la variance expliquée par le modèle

---
# Déviance unitaire, totale et nulle

La déviance unitaire est une mesure de la distance entre `$y$` et `$μ$`.

`$${\displaystyle d(y,y)=0}$$`
`$${\displaystyle d(y,\mu )>0\quad \forall y\neq \mu }$$`
???
Nous allons utiliser les déviances pour calculer pseudo R2.
--

La déviance totale `${\displaystyle D(\mathbf {y} ,{\hat {\boldsymbol {\mu }}})}$` d'un modèle avec prédictions `${\hat {\boldsymbol {\mu }}}$` de l'observation `$\mathbf {y}$` est la somme de ses déviances unitaires:

`$${\displaystyle D(\mathbf {y},{\hat {\boldsymbol {\mu }}})=\sum _{i}d(y_{i},{\hat {\mu }}_{i})}$$`

Maintenant, la déviance d'un modèle avec des prédictions `${\hat {\mu }}=E[Y|{\hat {\theta }}_{0}]$` peut être définie par sa **vraisemblance**:

`$$D(y,{\hat {\mu }})=2{\Big (}\log {\big (}p(y\mid {\hat {\theta }}_{s}){\big )}-\log {\big (}p(y\mid {\hat {\theta }}_{0}){\big )}{\Big )}$$`
avec les paramètres ajustées dans le modèle réduit ( `$\hat \theta_0$` ) et saturé ( `$\hat \theta_s$` ).

---
# Déviance unitaire, totale et nulle

La **déviance résiduelle** est définie comme 2 fois le logarithme du ratio de vraisemblance du modèle complet par rapport au modèle réduit. (La fonction ci-dessous est exactement la même que ci-dessus !)

`$$D(y,{\hat {\mu }})=2{\Big (}\log {\big (}p(\text{modèle saturé}){\big )}-\log {\big (}p(\text{modèle réduit}){\big )}{\Big )}$$`

Et, la **déviation nulle** est définie comme 2 fois le logarithme du ratio de vraisemblance du modèle complet par rapport au modèle nul (où les prédicteurs sont fixés à 1).

`$$D(y,{\hat {\mu }})=2{\Big (}\log {\big (}p(\text{modèle saturé}){\big )}-\log {\big (}p(\text{modèle nul}){\big )}{\Big )}$$`

???

Dans certaines conditions, si le modèle proposé décrit les données presque aussi bien que le modèle saturé, alors asymptotiquement D ∼ χ
2 K-p, où K et p sont le nombre de paramètres dans les modèles saturé et proposé, respectivement. Si le modèle proposé est mauvais, D sera plus grand que prévu par la distribution K-p de χ2.

---
# Déviance totale et nulle

Voici comment on peut comparer la déviance du modèle (déviance résiduelle) à la déviance d'un modèle nul (déviance nulle) dans `R`.

Le **modèle nul** est un modèle sans variables explicatives.

```R
null.model <- glm(Response.variable ~ 1, family = binomial)
```

Le **modèle à déviance saturée (ou totale)** est ici un modèle avec toutes les variables explicatives.

```r
full.model <- glm(response.variable ~ ., family = binomial)
```

---
# Pouvoir prédictif et ajustement du modèle

Dans `R`, nous pouvons extraire les déviances résiduelles et nulles directement à partir de l'`objet glm`:

.small[

```r
objects(logit.reg)
#  [1] "aic"               "boundary"          "call"             
#  [4] "coefficients"      "contrasts"         "control"          
#  [7] "converged"         "data"              "deviance"         
# [10] "df.null"           "df.residual"       "effects"          
# [13] "family"            "fitted.values"     "formula"          
# [16] "iter"              "linear.predictors" "method"           
# [19] "model"             "null.deviance"     "offset"           
# [22] "prior.weights"     "qr"                "R"                
# [25] "rank"              "residuals"         "terms"            
# [28] "weights"           "xlevels"           "y"
```
]

En utilisant les valeurs d'écart disponibles à partir de l'objet `logit.reg`, quelle est la valeur du pseudo `$R^2$` ?

???

Demandez aux participants d'écrire un code à partir de ces informations sur les déviances pour calculer manuellement `$R^2$` et écrire leur réponse dans le chat.

---
# Pouvoir prédictif et ajustement du modèle

Dans `R`, nous pouvons extraire les déviances résiduelles et nulles directement à partir de l'`objet glm`:

.small[

```r
pseudoR2 <- (logit.reg$null.deviance - logit.reg$deviance) / logit.reg$null.deviance
pseudoR2
# [1] 0.4655937
```
]

.comment[Donc, le modèle explique 46.6% de la variabilité des données.]

---
# Pouvoir prédictif et ajustement du modèle

Un pseudo-R2 de McFadden ajusté, qui pénalise pour le nombre de prédicteurs, peut être calculé comme suit :
.pull-left[
`\begin{equation}
R^2_{adj} = 1 - \frac{logL(M) - K}{logL(M_{null})}
\end{equation}`
]

.pull-right[
où `$K$` correspond au nombre supplémentaire de prédicteurs par rapport au modèle nul.
]

La qualité d'ajustement des modèles de régression logistique peut être exprimée par des variantes de statistiques pseudo-R2, telles que les mesures de Maddala (1983) ou de Cragg et Uhler (1970).

Lorsqu'on parle de *régressions logistiques*, les valeurs faibles de `$R^2$` sont souvent courantes.

---
# Pouvoir prédictif et ajustement du modèle

La fonction `DescTools::PseudoR2()` calcule plusieurs pseudo-R2.

La fonction R `DescTools::PseudoR2()` permet de calculer plusieurs pseudo-R2. En spécifiant `which = all`, calculez toutes les statistiques en même temps.

```r
logit.reg <- glm(pa ~ WatrCont + Topo, 
 data = mites, family = binomial(link = "logit"))
DescTools::PseudoR2(logit.reg, which = "all")
# McFadden McFaddenAdj CoxSnell Nagelkerke AldrichNelson 
# 0.4655937 0.3998373 0.4549662 0.6245898 0.3776866 
# VeallZimmermann Efron McKelveyZavoina Tjur AIC 
# 0.6674318 0.5024101 0.7064093 0.5114661 54.7623962 
# BIC logLik logLik0 G2 
# 61.5078819 -24.3811981 -45.6229593 42.4835224
```

???
Pour évaluer l'ajustement du modèle, les graphiques de diagnostique ne sont pas utiles. Il vaut mieux utiliser le test de Hosmer-Lemeshow qui évalue si les fréquences des événements observés correspondent ou non aux fréquences des événements attendus.

Vous pouvez le montrer aux étudiants avec 
library(vcdExtra)
HLtest(logit.reg)
---
exclude: true

### Exercice 2: Faire un test d'ajustement du modèle Hosmer-Lemeshow

Pour évaluer l'ajustement du modèle, les graphiques de diagnostique ne sont pas utiles. Il vaut mieux utiliser le [test de Hosmer-Lemeshow](https://en.wikipedia.org/wiki/Hosmer%E2%80%93Lemeshow_test):

Il s'agit en général de:
1. La division des données en groupes de taille similaire;
2. Ordonner sur la probabilité prédite (ou de manière équivalente, le prédicteur linéaire);
3. Comparer le nombre de réponses positives observées au nombre de réponses positives attendues dans chaque groupe;
4. Réalisation d'un test du chi carré sur le résultat.

Ce test évalue si les fréquences des événements observés correspondent ou non aux fréquences des événements attendus dans les sous-groupes de la population modèle.

.pull-left[

```r
library(vcdExtra)
HLtest(logit.reg)
# Hosmer and Lemeshow Goodness-of-Fit Test 
# 
# Call:
# glm(formula = pa ~ WatrCont + Topo, family = binomial(link = "logit"), 
#     data = mites)
#  ChiSquare df   P_value
#   3.421693  8 0.9051814
```
]
.pull-right[

```r
plot(HLtest(logit.reg))
```

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

.comment[Une valeur non significative indique un ajustement adéquat!]

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

1\. En utilisant le modèle créé avec le jeu de données `bacteria`, évaluez le pouvoir prédictif et l'ajustement de ce modèle.

2\. Comment améliorer le pouvoir explicatif du modèle?

???

Faites cet exercise en salles de réunion avec des groupes de 3-4 pour un durée de 5 minutes.

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

1\. En utilisant le modèle créé avec le jeu de données `bacteria`, évaluez le pouvoir prédictif et l'ajustement de ce modèle.

```r
null.d <- model.bact2$null.deviance
resid.d <- model.bact2$deviance
bact.pseudoR2 <- (null.d - resid.d) / null.d
bact.pseudoR2
# [1] 0.0624257
DescTools::PseudoR2(model.bact2, which = "all")
# McFadden McFaddenAdj CoxSnell Nagelkerke AldrichNelson 
# 0.06242570 0.02562309 0.05981734 0.09529510 0.05809758 
# VeallZimmermann Efron McKelveyZavoina Tjur AIC 
# 0.11689650 0.06151653 0.11098468 0.06275058 211.80606240 
# BIC logLik logLik0 G2 
# 225.38057258 -101.90303120 -108.68795256 13.56984272
```
--

2\. Comment améliorer le pouvoir explicatif du modèle?

*Ajouter des variables explicatives pertinentes pourrait certainement augmenter le pouvoir explicatif du modèle.*

*Mais, n'ayez pas peur de résultats non-significatifs!*

---
# GLM et données de proportions

Parfois, les données de proportions sont plus similaires à une régression logistique que ce que vous pensez...

En comptes discrets, nous pouvons, par exemple, mesurer le nombre de présences d'individus par rapport au nombre total de populations échantillonnées.

Nous obtiendrons ainsi un nombre proportionnel de "succès" dans l'observation des individus en divisant les comptes par les comptes totaux.

> Dans `glm()`, nous devons fournir des *poids a priori* si la variable de réponse est la proportion de succès.

---
# GLM et données de proportions

Les proportions peuvent être codés en fournissant le nombre de succès et des poids a priori dans la fonction:

```r
prop.reg <- glm(cbind(Galumna, totalabund - Galumna) ~ Topo + WatrCont,
 data = mites,
 family = binomial)
```

Les poids peuvent aussi être spécifiés dans `glm()`:

```r
prop.reg2 <- glm(prop ~ Topo + WatrCont,
 data = mites,
 family = binomial,
 weights = totalabund)
```
---
exclude: true

# Exercice 3

Dans R, on doit spécifier le nombre de fois qu'un événement s'est produit et le nombre de fois qu'un événement ne s'est pas produit:

```r
prop.reg <- glm(cbind(Galumna, totalabund - Galumna) ~ Topo + WatrCont, 
 data = mites, 
 family = binomial)
```

```r
summary(prop.reg)
```

---
exclude: true

# Exercice 3

.small[

```r
summary(prop.reg)
# 
# Call:
# glm(formula = cbind(Galumna, totalabund - Galumna) ~ Topo + WatrCont, 
#     family = binomial, data = mites)
# 
# Deviance Residuals: 
#     Min       1Q   Median       3Q      Max  
# -1.4808  -0.9699  -0.6327  -0.1798   4.1688  
# 
# Coefficients:
#              Estimate Std. Error z value Pr(>|z|)    
# (Intercept) -3.288925   0.422109  -7.792 6.61e-15 ***
# TopoHummock  0.578332   0.274928   2.104   0.0354 *  
# WatrCont    -0.005886   0.001086  -5.420 5.97e-08 ***
# ---
# Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# 
# (Dispersion parameter for binomial family taken to be 1)
# 
#     Null deviance: 140.702  on 69  degrees of freedom
# Residual deviance:  85.905  on 67  degrees of freedom
# AIC: 158.66
# 
# Number of Fisher Scoring iterations: 5
```
]

---
exclude: true

# Exercice 3

On peut coder le modèle directement avec les proportions:

```r
prop.reg2 <- glm(prop ~ Topo + WatrCont, data = mites,
 family = binomial, weights = totalabund)
```

---
class: inverse, center, middle

# GLM avec des données d'abondance

---
# Modéliser des données d'abondance

.large[Que sont des données d'abondance?]

Importez le jeu de données `faramea.csv` dans `R`.

```r
faramea <- read.csv('faramea.csv', header = TRUE)
```

Le nombre d'arbres de l'espèce *Faramea occidentalis* a été compté dans 43 quadrats sur l'île de Barro Colorado (Panama). Des données environnementales, comme l'élévation et la précipitation ont aussi été mesurées.

Examinons maintenant à quoi ressemble la distribution du nombre d'arbres par transect.

---
# Modéliser des données d'abondance

.large[Que sont des données d'abondance?]

Les données d'abondance sont charactérisées par:

- des valeurs positives : on ne peut pas compter **-7** individus
- des valeurs entières : on ne peut pas compter **7.56** individus
- une plus grande variance pour les grandes valeurs

<img src="workshop06-pres-fr_files/figure-html/unnamed-chunk-65-1.png" width="576" style="display: block; margin: auto;" />
???
Vous pouvez dire qu'en *moyenne*, vous avez 7.56 individus pour une situation spécifique. Mais la mesure est tout de même *discrète*.
---
# Modéliser des données d'abondance

.large[Comment modéliser des données d'abondance?]

L'élévation influence-t-elle l'abondance de *F. occidentalis*?

---
# Modéliser des données d'abondance

L'élévation influence-t-elle l'abondance de *F. occidentalis*?

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

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

Quelle distribution modèle le mieux ces observations?
1. Bernoulli
2. Normale
3. Poisson
4. Binomial

???

Ceci est un question de sondage à un choix. La réponse est 3. Poisson.

---
# Modéliser des données d'abondance

L'élévation influence-t-elle l'abondance de *F. occidentalis*?
.pull-left[
<img src="workshop06-pres-fr_files/figure-html/unnamed-chunk-69-1.png" width="432" style="display: block; margin: auto;" />
]

.pull-right[
<img src="workshop06-pres-fr_files/figure-html/unnamed-chunk-70-1.png" width="432" style="display: block; margin: auto;" />
]
.small[Quelle distribution modèle le mieux ces observations?
1. Bernoulli 
2. Normal 
3. **Poisson**
4. Binomial
]

---
# La distribution de Poisson

La distribution de poisson, qui spécifie la probabilité d'une variable aléatoire discrète Y, est données par :

`$$f(y, \,\mu)\, =\, Pr(Y = y)\, =\, \frac{\mu^y \times e^{-\mu}}{y!}$$`

`$$E(Y)\, =\, Var(Y)\, =\, \mu$$`

**Propriétés** :

- `$\mu$` est le paramètre de la distribution de Poisson
- spécifie la probabilité pour des valeurs entières uniquement
- la probabilité pour des valeurs négatives est nulle ( `$P(Y<0) = 0$`)
- moyenne = variance (permet l'hétérogénéité)

---
# Que se cache-t-il derrière un GLM Poisson?

Un GLM Poisson va modéliser la valeur de `$\mu$` en fonction de différentes variables explicatives.

.center[**Trois étapes**]

**Étape 1.** On suppose que `$Y_i$` suit une distribution de Poisson de moyenne et variance `$\mu_i$`

`$$Y_i = Poisson(\mu_i)$$`

`$$E(Y_i) = Var(Y_i) = \mu_i$$`

`$$f(y_i, \, \mu_i) = \frac{\mu^{y_i}_i \times e^{-\mu_i}}{y!}$$`

`$\mu_i$` correspond au nombre attendu d'individus

---
# Que se cache-t-il derrière un GLM Poisson?

**Étape 2.** On spécifie le prédicteur linéaire comme dans un modèle linéaire.

`$$\underbrace{\alpha}_\text{One intercept} + \underbrace{\beta}_\text{slope of 'Elevation'} \times \text{Elevation}_i$$`
--
**Étape 3.** La fonction de lien entre la moyenne de `$Y_i$` et la partie systématique est une fonction logarithmique et est écrite comme suit:

.pull-left[
 
`$$log(\mu_i) = \alpha + \beta \times \text{Elevation}_i$$`
.center[ou]

`$$\mu_i = e^{ \alpha + \beta \times \text{Elevation}_i}$$`
]

.small[
Ceci montre que l'impact de chaque variable explicative est multiplicatif.
Augmenter l'élévation de un augmente μ par le facteur exp( `$\beta_\text{Elevation}$` ).

Si `$βj = 0$`, alors `$exp(βj) = 1$` et `$μ$` n'est pas lié à `$x_j$`. Si `$βj > 0$` alors `$μ$` augmente si `$x_j$` augmente ; si `$βj < 0$` alors `$μ$` diminue si `$x_j$` augmente.
]

---
# Ajuster un GLM Poisson dans R

La fonction `glm()` permet de spécifier un GLM Poisson.

```r
glm.poisson = glm(Faramea.occidentalis ~ Elevation,
                  data = faramea, family = poisson)
```

L'argument `family` spécifie le type de distribution et la fonction lien.

Tout comme avec `lm()`, on peut accéder au résumé du modèle à l'aide de la fonction `summary()`

```r
summary(glm.poisson)
```

---
# Résumé du modèle

.pull-left2[
.small[

```r
summary(glm.poisson)
# 
# Call:
# glm(formula = Faramea.occidentalis ~ Elevation, family = poisson, 
# data = faramea)
# 
# Deviance Residuals: 
# Min 1Q Median 3Q Max 
# -3.3319 -2.7509 -1.5451 0.1139 11.3995 
# 
# Coefficients:
# Estimate Std. Error z value Pr(>|z|) 
# (Intercept) 1.7687001 0.1099136 16.092 < 2e-16 ***
# Elevation -0.0027375 0.0006436 -4.253 2.11e-05 ***
# ---
# Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# 
# (Dispersion parameter for poisson family taken to be 1)
# 
# Null deviance: 414.81 on 42 degrees of freedom
# Residual deviance: 388.12 on 41 degrees of freedom
# AIC: 462.01
# 
# Number of Fisher Scoring iterations: 10
```
]]

.pull-right2[
Estimés :

Intercept = `$\alpha$`

Élevation = `$\beta$`

]

.pull-right2[
 
Qu'en est-il de `Null deviance` et `Residual deviance`?!
]

???

Mentionner que la déviance résiduelle de 388.12 est beaucoup plus élevée que les degrés de liberté résiduels de 41 ce qui indique le surdispersion

---
exclude: false

# Estimation des paramètres

Dans notre modèle, les paramètres à estimer sont l'ordonnée à l'origine ( `$\alpha$` ) et le coefficient de régression de l'élevation ( `$\beta$` )

`$$log(\mu_i) = 1.769 - 0.0027 \times \text{Élevation}_i$$`

.center[ou]

`$$\mu_i = e^{1.769 - 0.0027 \times \text{Élevation}_i}$$`

---
exclude: false

# La déviance

Rappelez vous que pour estimer les paramètres inconnus, l'estimation par maximum de vraisemblance est utilisée.

La déviance résiduelle est approximativement la différence entre la vraisemblance d'un modèle saturé ($n$ paramètres pour chaque observation) et le modèle complet ($p$ paramètres):

`$$\text{Res dev} = 2 \, log(L(y;\,y)) - 2 \, log(L(y;\, \mu))$$`

Dans un GLM Poisson, la déviance résiduelle doit être égale au nombre de degrés de liberté résiduels.

.center[.alert[388.12 >> 41]]

> Notre déviance résiduelle est beaucoup plus grande que les degrés de liberté de notre modèle!

---
# La surdispersion

Pour une distribution Poisson, `$var[y] = μ$`. En pratique, cependant, 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 le
les événements qui sont comptés sont corrélés positivement.

.center[*Pouvez-vous penser à une situation en biologie qui peut provoquer la surdispersion?*]

--
 
**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 justifierait!

---
# La surdispersion

Quand la déviance résiduelle est supérieure au nombre de degrés de liberté résiduels, le modèle est **surdispersé**.

.pull-left[
On peut calculer le paramètre de surdispersion ( `$\phi$` ) par:
`$$\phi ~ = ~\frac{\text{Déviance résiduelle}}{\text{Degrés de liberté résiduels}}$$`
]

.pull-right[
Si, `$\phi$` est plus que `$1$`, alors nous devrons l'estimer avec:

`$$\widehat{\phi}=\frac{1}{n-k}\sum\frac{(Y_{i}-\hat{\mu_{i}})^2}{\hat{\mu_{i}}}$$`
]

```
# En R, vous pouvez obtenir phi avec:
sum(residuals(ton_obj_glm, type="pearson")^2)/df.residual(ton_obj_glm)
```

Puisque le GLM Poisson requiert que `$var[y] = μ$`, nous devons trouver une solution alternative.

.center[.large[**Solutions**]]

.pull-left[
1: Corriger la surdispersion en utilisant un **GLM quasi-Poisson**
]

.pull-right[
2: Choisir une autre distribution, comme **la negative binomiale**
]

???

where n is the sample size, k is the number of estimated parameters (including the intercept) and μiˆ=g−1(ηiˆ) is the fitted expectation of Yi (g−1 is the inverse link function). Note that the formulation above is simply the Pearson χ2 divided by the residual degrees of freedom: ϕˆ=χ2/df.

---
# GLM quasi-Poisson

La variance du modèle tient compte de la **surdispersion** en ajoutant le paramètre de surdispersion:

`$$E(Y_i) = \mu_i$$`

`$$Var(Y_I) = \phi \times \mu_i$$`
où `$\phi$` (*phi*) est le paramètre de dispersion. Il sera estimé avant les paramètres.

Le **prédicteur linéaire** et la **fonction de lien** restent les mêmes.

Corriger pour la surdispersion ne va pas affecter l'estimation des paramètres, mais leur **significativité**. En fait, les écarts-types des paramètres seront multipliés par `$\sqrt{\phi}$`.

.alert[Certaines p-values marginalement significatives peuvent devenir non significatives!]

---
exclude: false

# Ajuster un GLM quasi-Poisson dans `R`

Créez un nouveau GLM à l'aide de la famille `quasipoisson` ou mettez le modèle précédent à jour:

```r
glm.quasipoisson = glm(Faramea.occidentalis ~ Elevation, 
                       data = faramea,
                       family=quasipoisson)

glm.quasipoisson = update(glm.poisson, 
                          family = quasipoisson)
```

---
# Ajuster un GLM quasi-Poisson dans `R`

.pull-left2[
.small[

```r
summary(glm.quasipoisson)
# 
# Call:
# glm(formula = Faramea.occidentalis ~ Elevation, family = quasipoisson, 
#     data = faramea)
# 
# Deviance Residuals: 
#     Min       1Q   Median       3Q      Max  
# -3.3319  -2.7509  -1.5451   0.1139  11.3995  
# 
# Coefficients:
#              Estimate Std. Error t value Pr(>|t|)    
# (Intercept)  1.768700   0.439233   4.027 0.000238 ***
# Elevation   -0.002738   0.002572  -1.064 0.293391    
# ---
# Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# 
# (Dispersion parameter for quasipoisson family taken to be 15.96936)
# 
#     Null deviance: 414.81  on 42  degrees of freedom
# Residual deviance: 388.12  on 41  degrees of freedom
# AIC: NA
# 
# Number of Fisher Scoring iterations: 10
```
]]
.pull-right2[
**Mêmes estimés mais**

.small[Les écarts-types des paramètres sont multipliés par]

`$$\sqrt{\phi} = 4$$`

`0.0006436 * 4 = 0.00257`

<- `$\phi$`

<- .small[L'AIC n'est pas défini!]
]

---
# Ajuster un GLM quasi-Poisson sous R

Testons l'effet de l'élévation par une analyse de déviance.

```r
null.model <- glm(Faramea.occidentalis ~ 1, 
 data = faramea,
 family = quasipoisson)
anova(null.model, glm.quasipoisson, test = "Chisq")
# Analysis of Deviance Table
# 
# Model 1: Faramea.occidentalis ~ 1
# Model 2: Faramea.occidentalis ~ Elevation
# Resid. Df Resid. Dev Df Deviance Pr(>Chi)
# 1 42 414.81 
# 2 41 388.12 1 26.686 0.1961
```

---
# Paramètre de dispersion

.center[![:scale 80%](images/dispParam.png)]

???

Un GLM quasi-Poisson peut corriger un peu de surdispersion, mais une autre distribution pourrait être approprié si il y en a trop.

---
# GLM binomiale négative

Une distribution binomiale négative est favorable quand la surdispersion est forte.

- La distribution a **deux paramètres** `$\mu$` and `$k$`. `$k$` contrôle pour la dispersion (plus la dispersion est forte, plus `$k$` est petit);
- C'est une combinaison de deux distributions (**Poisson** et **gamma**);
- Les `$Y_i$` suivent une distribution de Poisson dont la moyenne `$\mu$` suit une distribution Gamma!

Les valeurs prédites suivent:

`$$E(Y_i) = \mu_i$$`

et la fonction de variance:

`$$Var(Y_i) = \mu_i + \frac{\mu^2_i}{k}$$`

---
# Ajuster une binomiale négative dans `R`

NB La distribution binomiale n'est pas dans la fonction `glm()`, donc il faut installer et charger la paquet `MASS`.

```r
install.packages('MASS')
```

```r
glm.negbin <- glm.nb(Faramea.occidentalis ~ Elevation, data = faramea)
```

```r
summary(glm.negbin)
```

---
# Ajuster une binomiale négative dans `R`

.pull-left2[
.small[

```r
summary(glm.negbin)
# 
# Call:
# glm.nb(formula = Faramea.occidentalis ~ Elevation, data = faramea, 
#     init.theta = 0.2593107955, link = log)
# 
# Deviance Residuals: 
#      Min        1Q    Median        3Q       Max  
# -1.36748  -1.17564  -0.51338  -0.05226   2.25716  
# 
# Coefficients:
#              Estimate Std. Error z value Pr(>|z|)    
# (Intercept)  2.369226   0.473841    5.00 5.73e-07 ***
# Elevation   -0.007038   0.002496   -2.82  0.00481 ** 
# ---
# Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# 
# (Dispersion parameter for Negative Binomial(0.2593) family taken to be 1)
# 
#     Null deviance: 41.974  on 42  degrees of freedom
# Residual deviance: 36.343  on 41  degrees of freedom
# AIC: 182.51
# 
# Number of Fisher Scoring iterations: 1
# 
# 
#               Theta:  0.2593 
#           Std. Err.:  0.0755 
# 
#  2 x log-likelihood:  -176.5090
```
]]
.pull-right2[

`theta` `$= k$`
]

---
# Visualiser le modèle final

**Étape 1** Représenter les données et utiliser les estimations des paramètres pour visualiser le modèle.

`$$\mu_i = e^{2.369 - 0.007 \times Elevation_i}$$`

Utilisez `summary()` pour obtenir les paramètres

```r
summary(glm.negbin)$coefficients[1, 1]
summary(glm.negbin)$coefficients[2, 1]
```

---
# Visualiser le modèle final

**Étape 2** Utilisez les écarts-types pour construire l'intervalle de confiance.

```r
summary(glm.negbin)$coefficients[1, 2]
summary(glm.negbin)$coefficients[2, 2]
```

`$$\text{Limite sup} = e^{[\alpha - 1.96 \times SE_{\alpha}] + [\beta - 1.96 \times SE_{\beta}] \times \text{Élevation}_i}$$`

`$$\text{Limit inf} = e^{[\alpha + 1.96 \times SE_{\alpha}] + [\beta + 1.96 \times SE_{\beta}] \times \text{Élevation}_i}$$`

---
# Visualiser le modèle final

```r
pp <- predict(glm.negbin, 
 newdata = data.frame(Elevation = 1:800), 
 se.fit = TRUE)
linkinv <- family(glm.negbin)$linkinv ## fonction lien-inverse
pframe$pred0 <- pp$fit
pframe$pred <- linkinv(pp$fit)
sc <- abs(qnorm((1-0.95)/2)) ## Normal approx. to likelihood
pframe <- transform(pframe, 
 lwr = linkinv(pred0-sc*pp$se.fit), 
 upr = linkinv(pred0+sc*pp$se.fit))
# sinon, utiiser predict() avec type="response"

plot(faramea$Elevation, faramea$Faramea.occidentalis, 
     ylab = 'Number of F. occidentalis', xlab = 'Elevation(m)')
lines(pframe$pred, lwd = 2)
lines(pframe$upr, col = 2, lty = 3, lwd = 2)
lines(pframe$lwr, col = 2, lty = 3, lwd = 2)
```

---
# Visualiser le modèle final

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

Utilisez le jeu de données `mites`! Modélisez l'abondance de l'espèce *Galumna* en fonction des caractéristiques du substrat (son contenu en eau `WatrCont` et sa densité `SubsDens`).

- Faut-il contrôler pour la surdispersion?
- Quelles variables explicatives ont un effet significatif?
- Selectionnez le meilleur modèle!

```r
mites <- read.csv("mites.csv", header = TRUE)
```

???

Faites cet exercice en salles de réunion avec des equipes de 3-4, pour une durée de 10 minutes. N'oubliez pas de montrer les conseils sur la prochaine diapositive.

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

Retirez une variable à la fois et comparez le modèle imbriqué au modèle saturé (ou complet):

```r
drop1(MyGLM, test = "Chi")
```

Spécifiez un modèle imbriqué manuellement, appelez le `MyGLM2`, et utilisez la fonction `anova()`:

```r
anova(MyGLM, MyGLM2, test = "Chi")
```

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

.small[

```r
# GLM Poisson
glm.p = glm(Galumna~WatrCont+SubsDens, data=mites, family=poisson)
# GLM quasi-Poisson
glm.qp = update(glm.p,family=quasipoisson)
# sélection du modèle
drop1(glm.qp, test = "Chi")
# Single term deletions
# 
# Model:
# Galumna ~ WatrCont + SubsDens
# Df Deviance scaled dev. Pr(>Chi) 
# <none> 101.49 
# WatrCont 1 168.10 31.711 1.789e-08 ***
# SubsDens 1 108.05 3.125 0.07708 . 
# ---
# Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
```
```r
# ou
glm.qp2 = glm(Galumna~WatrCont, data=mites, family=quasipoisson)
anova(glm.qp2, glm.qp, test="Chisq")
```
]

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

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

---
# Autres distributions

- **Transformation logit des données** souvent utilisée avec `lm()` pour les pourcentages et les proportions quand la distribution binomiale n'est pas appropriée. Quand non selectionné à partir de quantités fixées (e.g. pourcentage de couverture, grades scolaires, etc).
- **Distribution log-normal dans un glm**, évite d'avoir à log-transformer les données.
- **Distribution Gamma**. Similaire à une log-normal, plus flexible.
- **Distribution tweedie**. Famille de distributions flexible. Utile pour des données avec un mélange de 0 et de valeurs positives (pas forcément des comptes).
- **Poisson ou negative binomiale à inflation de zéro**. Quand les données comprennent un nombre excessif de zéros, venant d'un processus différent de celui qui génère les comptes.

---

# Ressources additionnelles sur les GLMs

**Livres** :

- B. Bolker (2009) Ecological Models and Data in R. Princeton University Press.
- A. Zuur et al. (2009) Mixed Effects Models and Extensions in Ecology with R. Springer.

**Articles** :

- [Harrison et al. (2018), PeerJ, DOI 10.7717/peerj.4794
](http://dx.doi.org/10.7717/peerj.4794)

**Sites internet** :

- GLMM for ecologists (http://glmm.wikidot.com) .small[.comment[A great website on GLMM with a Q&A section!]]

---
class: inverse, center, bottom

# Merci pour votre participation à cet atelier!

<hr>

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