Atelier 8: Modèles additifs généralisés

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

---

# À propos de cet atelier

[![badge](https://img.shields.io/static/v1?style=for-the-badge&label=Présentation&message=08&color=BF616A)](https://r.qcbs.ca/workshop08/pres-fr/workshop08-pres-fr.html)    [![badge](https://img.shields.io/static/v1?style=for-the-badge&label=Script&message=08&color=D08770&logo=r)](https://r.qcbs.ca/workshop08/book-fr/workshop08-script-fr.R) 
[![badge](https://img.shields.io/static/v1?style=for-the-badge&label=Livre&message=08&color=EBCB8B)](https://r.qcbs.ca/workshop08/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-08/)
[![badge](https://img.shields.io/static/v1?style=for-the-badge&label=GitHub&message=08&color=B48EAD&logo=github)](https://github.com/QCBSRworkshops/workshop08)

---

# Librairies et jeux de données

Pour cet atelier, nous travaillerons avec les jeux de données suivants :

* [ISIT.csv](https://r.qcbs.ca/workshop08/pres-en/data/ISIT.csv)

Vous devriez également vous assurer que vous avez téléchargé, installé et chargé les librairies R suivants:

* [ggplot2](https://cran.r-project.org/package=ggplot2)
* [itsadug](https://cran.r-project.org/package=itsadug)
* [mgcv](https://cran.r-project.org/package=mgcv)

```R
install.packages(c('ggplot2', 'itsadug', 'mgcv'))
```

---

# Prérequis

Pour suivre cet atelier, nous recommandons:

+ Assez d'expérience avec le logiciel R pour exécuter un script et examiner les données et les objets R;
+ Une connaissance de base de la régression linéaire.

---
# Aperçu de l'atelier

Nous commencerons par les bases de comment spécifier et implémenter des GAMs en R:

1. Le modèle linéaire... et où il échoue
2. Introduction aux GAMs
3. Le fonctionnement des GAMs
4. GAM avec plusieurs termes non-linéaires
5. GAM avec des termes d'interaction

Nous discuterons ensuite de la généralisation du modèle additif, incluant:

6. La validation d'un GAM
7. Choisir une autre distribution
8. Changer la fonction de base
9. Introduction rapide aux GAMMs

---
# Objectifs d'apprentissage

1. Utiliser la librairie `mgcv` pour modéliser les relations non linéaires,
2. Évaluer la sortie d'un Modèle Additif Généralisé (GAM) afin de mieux comprendre nos données,
3. Utiliser des tests pour déterminer si nos relations correspondent à des modèles non linéaires ou linéaires,
4. Ajouter des interactions non linéaires entre les variables explicatives,
5. Comprendre l'idée d'une fonction de base (basis function) et pourquoi ça rend les GAMs si puissants,
6. Comment modéliser la dépendance dans les données (autocorrélation, structure hiérarchique) en utilisant les GAMMs.

---
class: inverse, center, middle
# 1. Le modèle linéaire

<hr>

## ... et où il échoue

---
# La régression linéaire

La régression linéaire est parmi les méthodes les plus performantes et plus utilisées. Elle nous permet de modéliser une variable réponse en fonction de facteurs prédictifs et d'une erreur résiduelle.

Comme on a vu dans l'[Atelier 4: Modèles linéaires](https://r.qcbs.ca/fr/workshops/r-workshop-04/), le modèle linéaire fait quatre suppositions importantes:

1. Relation linéaire entre les variables de réponse et les variables prédicteurs:  `$$y_i = \beta_0 + \beta_1 \times x_i + \epsilon_i$$`
2. L'erreur est distribuée normalement: `$$\epsilon_i \sim \mathcal{N}(0,\,\sigma^2)$$`
3. La variance des erreurs est constante
4. Chaque erreur est indépendante des autres (homoscédasticité)

*Modèle linéaire avec plusieurs prédicteurs:*

`$$y_i = \beta_0 + \beta_1x_{1,i}+\beta_2x_{2,i}+\beta_3x_{3,i}+...+\beta_kx_{k,i} + \epsilon_i$$`
???
Un modèle linéaire peut parfois s'adapter à certains types de réponses non linéaires (par exemple `$x^2$`), mais cette approche repose fortement sur des décisions qui peuvent être soit arbitraires, soit bien informées, et est beaucoup moins flexible que l'utilisation d'un modèle additif.

---
# La régression linéaire

Les modèles linéaires fonctionnent très bien dans certains cas spécifiques où tous ces critères sont respectés:

---
# La régression linéaire

En réalité, il est souvent impossible de respecter ces critères. Dans de nombreux cas, les modèles linéaires sont inappropriés:

---
# La régression linéaire

**Quel est le problème et comment le régler?**

Un **modèle linéaire** essaye d'ajuster la meilleure **droite** qui passe au milieu des données, cela ne fonctionne donc pas pour tous les jeux de données.

En revanche, les **GAM** font cela en ajustant une **fonction de lissage non-linéaire** à travers les données, mais tout en contrôlant le degré de courbure de la ligne (*plus d'informations suivront*).

---
class: inverse, center, middle

## 2. Introduction aux GAMs
---
# Modèle Additif Généralisé (GAM)

Utilisons un exemple pour démontrer la différence entre une régression linéaire et un modèle additif. Nous allons utiliser le jeu de données `ISIT`.

```r
isit <- read.csv("data/ISIT.csv")
head(isit)
```

> Ce jeu de donnée comporte des mesures de bioluminescence en relation à la profondeur (*depth*), la station de rechercher et la saison (*Season*).

Prenons que les données de la deuxième saison pour l'instant:

```r
isit2 <- subset(isit, Season==2)
```

---
# Modèle Additif Généralisé (GAM)

Commençons par essayer d'ajuster un modèle de régression linéaire à la relation entre `Sources` et `SampleDepth`.

```r
linear_model <- gam(Sources ~ SampleDepth, data = isit2)
summary(linear_model)
...
# 
# Parametric coefficients:
# Estimate Std. Error t value Pr(>|t|) 
# (Intercept) 30.9021874 0.7963891 38.80 <2e-16 ***
# SampleDepth -0.0083450 0.0003283 -25.42 <2e-16 ***
# 
# 
# R-sq.(adj) = 0.588 Deviance explained = 58.9%
# GCV = 60.19 Scale est. = 59.924 n = 453
...
```

Le modèle linéaire explique une bonne partie de la variance de notre jeu de données ( `$R_{adj}$` = 0.588), ce qui veut dire que notre modèle est super bon, non?

???

Nous pouvons utiliser la commande `gam()` de la librairie `mgcv` pour modéliser une régression par les moindres carrés. Nous verrons plus loin comment utiliser `gam()` pour spécifier un terme lissé et non linéaire.

---
# Modèle Additif Généralisé (GAM)

Voyons comment notre modèle cadre avec les données:

```r
data_plot <- ggplot(data = isit2, aes(y = Sources, x = SampleDepth)) + 
 geom_point() +
 geom_line(aes(y = fitted(linear_model)),
 colour = "red", size = 1.2) + 
 theme_bw()
data_plot
```

???

Les suppositions de la régression linéaire sont-elles satisfaites dans ce cas? Comme vous l'avez peut-être remarqué, nous ne respectons pas les conditions du modèle linéaire:

1.  Il existe une forte relation _non linéaire_ entre `Sources` et `SampleDepth`.
2.  L'erreur n'est _pas_ normalement distribuée.
3.  La variance de l'erreur n'est _pas_ homoscédastique.
4.  Les erreurs ne sont _pas_ indépendantes les unes des autres.

---
exclude:true
# Modèles Additifs Généralisés (GAMs)

Examinons un exemple! Premièrement, nous allons générer des données et les représenter graphiquement.

```r
library(ggplot2)
set.seed(10)
n <- 250
x <- runif(n,0,5)
y_model <- 3*x/(1+2*x)
y_obs <- rnorm(n,y_model,0.1)
data_plot <- qplot(x, y_obs) +
 geom_line(aes(y=y_model)) +
 theme_bw()
data_plot
```

---
exclude:true
# GAM

Si nous modélisions cette relation par une régression linéaire, les résultats ne respecteraient pas les suppositions énumérées ci-dessus.

---
# Modèles Additifs Généralisés (GAMs)

**Relation entre la variable réponse et le prédicteur**

Une variable prédicteur:
`$$y_i = \beta_0 + f(x_i) + \epsilon$$`

Plusieurs variables prédicteurs:
`$$y_i = \beta_0 + f_1(x_{1,i}) + f_2(x_{2,i}) + ... + \epsilon$$`

Un des grands avantages d'utiliser un GAM est que la forme optimale de la non-linéarité, i.e. **le degré de lissage** de `$f(x)$` est contrôlée en utilisant une régression pénalisée qui est déterminée automatiquement est déterminée automatiquement selon la méthode d'ajustement (généralement le *maximum de vraisemblance* ou *maximum likelihood*).

???

Au sens strict, les équations concernent un GAM gaussien avec lien d'identité, qui est aussi appelé "modèle additif" (sans "généralisé").

Étant donné que la fonction de lissage `$f(x_i)$` est non linéaire et
locale, l'ampleur de l'effet de la variable explicative peut varier en
fonction de la relation entre la variable et la réponse.

Autrement dit, contrairement à un coefficient fixe `$\beta x_i$`, la fonction `$f$`
peut changer tout au long du gradient `$x_i$`. Le degré de lissage de `$f$`
est contrôlée en utilisant une régression pénalisée qui est déterminée
automatiquement à l'aide d'une méthode de validation croisée
généralisée (GCV) de la librairie `mgcv`.

---
# Modèles Additifs Généralisés (GAMs)

Essayons de modéliser les données à l'aide d'une fonction de lissage `s(x)` avec `mgcv::gam()`

```r
gam_model <- gam(Sources ~ s(SampleDepth), data = isit2)
```
--

```
...
# Family: gaussian 
# Link function: identity 
# 
# Parametric coefficients:
# Estimate Std. Error t value Pr(>|t|) 
# (Intercept) 12.8937 0.2471 52.17 <2e-16 ***
# 
# Approximate significance of smooth terms:
# edf Ref.df F p-value 
# s(SampleDepth) 8.908 8.998 214.1 <2e-16 ***
# ---
# Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# 
# R-sq.(adj) = 0.81 Deviance explained = 81.4%
# GCV = 28.287 Scale est. = 27.669 n = 453
...
```

???

Nous pouvons essayer de construire un modèle plus approprié en ajustant les données avec un terme lissé (non-linéaire). Dans `mgcv::gam()`, les termes lissés sont spécifiés par des expressions de la forme `s(x)`, où `$x$` est la variable prédictive non linéaire que nous voulons lisser. Dans ce cas, nous voulons appliquer une fonction de lissage à `SampleDepth`.

La variance expliquée par notre modèle a augmenté de plus de 20% ($R_{adj}$ = 0.81)!

---
# Modèles Additifs Généralisés (GAMs)

Note: contrairement à un coefficient fixe `$\beta$`, la fonction de lissage peut changer tout au long du gradient `$x$`.

???

Lorsque nous comparons l'ajustement des modèles linéaire (rouge) et non linéaire (bleu), il est clair que ce dernier cadre mieux avec nos données

---
exclude:true
# GAM

Essayons de modéliser les données à l'aide d'une fonction de lissage `s(x)` avec `mgcv::gam()`

```r
library(mgcv)
gam_model <- gam(y_obs ~ s(x))
summary(gam_model)

data_plot <- data_plot +
 geom_line(colour = "blue", size = 1.2, aes(y = fitted(gam_model)))
data_plot
```

---
exclude:true
# GAM

.comment[Note: contrairement à un coefficient fixe `$\beta$`, la fonction de lissage peut changer tout au long du gradient `$x$`]

---
# GAM

La librairie `mgcv` comprend également une fonction `plot` qui, par
défaut, nous permet de visualiser la non-linéarité du modèle.

```r
plot(gam_model)
```

---
# Test de linéarité avec GAM

Comment tester si le modèle non linéaire offre une amélioration significative par rapport au modèle linéaire?

On peut utiliser `gam()` et `AIC()` pour comparer la performance d'un modèle linéaire contenant `x` comme prédicteur linéaire à la performance d'un modèle non linéaire contenant `s(x)` comme prédicteur lisse.

```r
linear_model <- gam(Sources ~ SampleDepth, data = isit2)
smooth_model <- gam(Sources ~ s(SampleDepth), data = isit2)
AIC(linear_model, smooth_model)
# df AIC
# linear_model 3.00000 3143.720
# smooth_model 10.90825 2801.451
```

Ici, l'AIC du GAM lissé est plus bas, ce qui indique que l'ajout d'une fonction de lissage améliore la performance du modèle. _La linéarité n'est donc pas soutenue par nos données._

???

En d'autres termes, on demande si l'ajout d'une fonction lisse au modèle linéaire améliore l'ajustement du modèle à nos données.

Pour expliquer brièvement, le critère d'information d'Akaike (AIC) est une mesure comparative de la performance d'un modèle, où des valeurs plus basses indiquent qu'un modèle est "plus performant" par rapport aux autres modèles considérés.

---
exclude:true
# Test de linéarité avec GAM

On peut utiliser `gam()` et `AIC()` pour tester si une supposition de linéarité est justifiée. Ici, on demande si l'ajout d'une fonction lisse au modèle linéaire améliore l'ajustement du modèle à nos données:

```r
linear_model <- gam(y_obs ~ x) # ajuster un modèle linéaire régulier avec gam()
nested_gam_model <- gam(y_obs ~ s(x) + x)
AIC(linear_model, nested_gam_model, test = "Chisq")
```

---
exclude:true
# Test de linéarité avec GAM

```r
linear_model <- gam(y_obs ~ x) # ajuster un modèle linéaire régulier avec gam()
nested_gam_model <- gam(y_obs ~ s(x) + x)
anova(linear_model, nested_gam_model, test = "Chisq")
```

.comment[Notez que le modèle `y_obs ~ s(x)` donne exactement les même résultats que `y_obs ~ s(x) + x`. Nous utilisons `$s(x) + x$` pour illustrer l'imbrication du modèle, mais le `$+ x$` peut être omis.]

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

Essayons maintenant de déterminer si les données enregistrées lors de la première saison doivent être modélisées par une régression linéaire ou par un modèle additif.

Répétons le test de comparaison avec `gam()` et `AIC()` en utilisant les données de la première saison seulement:

```r
isit1 <- subset(isit, Season == 1)
```

1. Ajustez un modèle linéaire et un GAM à la relation entre `Sources` et `SampleDepth`.
2. Déterminez si l'hypothèse de linéarité est justifiée pour ces données.
3. Quels sont les degrés de liberté effectifs du terme non-linéaire?

???

Nous n'avons pas encore discuté des degrés de liberté effectifs (**EDF**), mais ils sont un outil clé pour nous aider à interpréter l'ajustement d'un GAM. Gardez ce terme en tête. Plus sur ce sujet dans les prochaines sections!

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

__1.__ Ajustez un modèle linéaire et un GAM à la relation entre `Sources` et `SampleDepth`.

```r
linear_model_s1 <- gam(Sources ~ SampleDepth, data = isit1)
smooth_model_s1 <- gam(Sources ~ s(SampleDepth), data = isit1)
```

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

__2.__ Déterminez si l'hypothèse de linéarité est justifiée pour ces données.

Comme ci-dessus, la visualisation de la courbe du modèle sur notre ensemble de données est une excellente première étape pour déterminer si notre modèle est bien conçu.

```r
ggplot(isit1, aes(x = SampleDepth, y = Sources)) +
  geom_point() +
  geom_line(colour = "red", size = 1.2,
            aes(y = fitted(linear_model_s1))) +
  geom_line(colour = "blue", size = 1.2,
            aes(y = fitted(smooth_model_s1))) +
  theme_bw()
```

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

__2.__ Déterminez si l'hypothèse de linéarité est justifiée pour ces données.

Comme ci-dessus, la visualisation de la courbe du modèle sur notre ensemble de données est une excellente première étape pour déterminer si notre modèle est bien conçu.

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

On peut compléter cela par une comparaison quantitative des performances du modèle en utilisant `AIC()`.

```r
AIC(linear_model_s1, smooth_model_s1)
#                       df      AIC
# linear_model_s1 3.000000 2324.905
# smooth_model_s1 9.644938 2121.249
```

> Le score AIC moins élevé indique que le modèle lissé est plus performant que le modèle linéaire, ce qui confirme que la linéarité n'est pas appropriée pour notre ensemble de données.

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

__3.__ Quels sont les degrés de liberté effectifs du terme non-linéaire?

Pour obtenir les degrés de liberté effectifs, il suffit d'imprimer notre objet du modèle:

```r
smooth_model_s1
# 
# Family: gaussian 
# Link function: identity 
# 
# Formula:
# Sources ~ s(SampleDepth)
# 
# Estimated degrees of freedom:
# 7.64  total = 8.64 
# 
# GCV score: 32.13946
```

Les degrés de liberté effectifs (EDF) sont >> 1.

> Gardez cela à l'esprit, car nous reviendrons bientôt sur les EDF!

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

Nous allons maintenant essayer cela avec d'autres données générées aléatoirement.

```r
n <- 250
x_test <- runif(n, -5, 5)
y_test_fit <- 4 * dnorm(x_test)
y_test_obs <- rnorm(n, y_test_fit, 0.2)
```

1. Ajustez un modèle linéaire et un GAM à la relation entre `x_test` et `y_test_obs`.
2. Déterminez si l'hypothèse de linéarité est justifiée pour ces données.
3. Quels sont les degrés de liberté effectifs du terme non-linéaire ?

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

```r
linear_model_test <- gam(y_test_obs ~ x_test)
nested_gam_model_test <- gam(y_test_obs ~ s(x_test) + x_test)

AIC(linear_model_test, nested_gam_model_test, test="Chisq")
```

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

```r
qplot(x_test, y_test_obs) +
  geom_line(aes(y = y_test_fit)) +
  theme_bw()
```

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

```r
nested_gam_model_test
```

**Réponse** Oui la non-linéarité est justifiée. Les degrés de liberté effectifs (**EDF**) sont >> 1 (on reviendra la dessus bientôt).

---
class: inverse, center, middle

## 3. Le fonctionnement des GAMs

---
# Le fonctionnement des GAMs

Nous allons maintenant prendre quelques minutes pour regarder comment fonctionnent les GAMs. Commençons en considérant d'abord un modèle qui contient une fonction lisse `$f$` d'une covariable, `$x$` :

`$$y_i = f(x_i) + \epsilon_i$$`

Pour estimer la fonction `$f$`, nous avons besoin de représenter l'équation ci-dessus de manière à ce qu'elle devienne un modèle linéaire. Cela peut être fait en définissant des fonctions de base, `$b_j(x)$`, dont est composée `$f$` :

`$$f(x) = \sum_{j=1}^q b_j(x) \times \beta_j$$`

---
# Exemple: une base polynomiale

Supposons que `$f$` est considérée comme un polynôme d'ordre 4, de sorte que l'espace des polynômes d'ordre 4 et moins contient `$f$`. Une base de cet espace serait alors :

`$$b_0(x)=1 \ , \quad b_1(x)=x \ , \quad b_2(x)=x^2 \ , \quad b_3(x)=x^3 \ , \quad b_4(x)=x^4$$`

Alors `$f(x)$` devient :

`$$f(x) = \beta_0 + x\beta_1 +  x^2\beta_2 + x^3\beta_3 + x^4\beta_4$$`

... et le modèle complet devient :

`$$y_i = \beta_0 + x_i\beta_1 +  x^2_i\beta_2 + x^3_i\beta_3 + x^4_i\beta_4 + \epsilon_i$$`

---
# Exemple: une base polynomiale

Chaque fonction de base est multipliée par un paramètre à valeur réelle, `$\beta_j$`, et est ensuite additionnée pour donner la courbe finale `$f(x)$`.

En faisant varier le coefficient `$\beta_j$`, on peut faire varier la forme de `$f(x)$` pour produire une fonction polynomiale d'ordre 4 ou moins.

---
# Exemple: une base de spline cubique

Un spline cubique est une courbe construite à partir de sections d'un polynôme cubique reliées entre elles de sorte qu'elles sont continues en valeur. Chaque section du spline a des coefficients différents.

---
# Exemple: une base de spline cubique

Voici une représentation d'une fonction lisse utilisant une base de régression spline cubique de rang 5 avec des nœuds situés à incréments de 0.2:

Dans cet exemple, les nœuds sont espacés uniformément à travers la gamme des valeurs observées de x. Le choix du degré de finesse du modèle est pré-déterminé par le nombre de nœuds, qui était arbitraire.

---
#### Contrôler le degré de lissage avec des splines de régression pénalisés

Au lieu de contrôler le lissage (non linéarité) en modifiant le nombre de nœuds, nous gardons celui-ci fixé à une taille un peu plus grande que raisonnablement nécessaire et on contrôle le lissage du modèle en ajoutant une pénalité sur le niveau de courbure. Donc, plutôt que d'ajuster le modèle en minimisant (comme avec la méthode des moindres carrés) :

`$$||y - XB||^{2}$$`

Le modèle peut être ajusté en minimisant:

`$$||y - XB||^{2} + \lambda \int_0^1[f^{''}(x)]^2dx$$`

Quand `$\lambda$` tend vers `$∞$`, le modèle devient linéaire.

---
#### Contrôler le degré de lissage avec des splines de régression pénalisés

Si `$\lambda$` est trop élevé, les données seront trop lissées et si elle est trop faible, les données ne seront pas assez lissées. Idéalement, il serait bon de choisir une valeur `$\lambda$` de sorte que le `$\hat{f}$` prédit est aussi proche que possible du `$f$` observé.

Un critère approprié pourrait être de choisir `$\lambda$` pour minimiser :

`$$M = 1/n \times \sum_{i=1}^n (\hat{f_i} - f_i)^2$$`

Étant donné que `$f$` est inconnue, `$M$` doit être estimé.

Les méthodes recommandées pour ce faire sont le maximum de vraisemblance (maximum likelihood, *ML*) ou l'estimation par maximum de vraisemblance restreint (restricted maximum likelihood, *REML*).

La validation croisée généralisée (*GCV*) est une autre possibilité.

---
exclude:true
# Principe de validation croisée

1. ajustement faible par rapport aux données et ne fait pas mieux avec le point manquant.

2. très bon ajustement de la courbe du signal sous-jacent, le lissage passe à travers le bruit et la donnée manquante est plutôt bien prédite.

3. la courbe ajuste le bruit aussi bien que le signal, la variabilité supplémentaire amène à prédire la donnée manquante plutôt mal.

---
exclude:true
# Principe de validation croisée

---

## 4. GAM avec plusieurs termes non-linéaires

---

# GAM à variables linéaires et non linéaires

Avec les GAMs, il est facile d'ajouter des termes non linéaires et linéaires dans un seul modèle, plusieurs termes non linéaires ou même des interactions non linéaires.

Dans cette section, nous allons utiliser les données de `ISIT` de nouveau.

Nous allons essayer de modéliser la réponse `Sources` avec les prédicteurs `Season` and `SampleDepth` simultanément.

Tout d'abord, nous devons convertir notre prédicteur qualitatif (`Season`) en facteur:

```r
head(isit)
isit$Season <- as.factor(isit$Season)
```

???

Rappelez-vous de ce jeu de données présenté dans les sections précédentes? Le jeu de données ISIT est composé des niveaux de bioluminescence (`Sources`) en fonction de la profondeur, des saisons et des différentes stations.

---
# GAM à variables linéaires et non linéaires

Commençons par un modèle de base comprenant un terme non linéaire (`SampleDepth`) et un facteur qualitatif (`Season` avec 2 niveaux).

```r
basic_model <- gam(Sources ~ Season + s(SampleDepth), data = isit, method = "REML")
basic_summary <- summary(basic_model)
```

Le tableau `p.table` donne des informations sur les termes linéaires:

```r
basic_summary$p.table
#             Estimate Std. Error  t value     Pr(>|t|)
# (Intercept) 7.253273  0.3612666 20.07734 1.430234e-72
# Season2     6.156130  0.4825491 12.75752 5.525673e-34
```

Le tableau `s.table` nous donne donne des informations sur le terme non linéaire:

```r
basic_summary$s.table
#                     edf   Ref.df        F p-value
# s(SampleDepth) 8.706426 8.975172 184.3583       0
```

---
exclude:true
# GAM à plusieurs variables

Avec les GAMs, il est facile d'ajouter des termes non linéaires et linéaires dans un seul modèle, plusieurs termes non linéaires ou même des interactions non linéaires.

Dans cette section, nous allons utiliser un ensemble de données générées automatiquement par `mgcv::gamSim()`.

```r
# ?gamSim
gam_data <- gamSim(eg = 5)
head(gam_data)
```

Nous allons essayer de modéliser la réponse `y` avec les prédicteurs `x0` à `x3`.

---
exclude:true
# GAM à plusieurs variables

Commençons par un modèle de base comprenant un terme non linéaire (`x1`) et un facteur qualitatif (`x0` avec 4 niveaux).

```r
basic_model <- gam(y ~ x0 + s(x1), data = gam_data)
basic_summary <- summary(basic_model)
basic_summary$p.table

basic_summary$s.table
```
.comment[La sortie de `p.table` fournit le tableau de résultats pour chaque terme paramétrique

Le tableau `s.table` nous donne les résultats du terme non linéaire.
]

---
# GAM à variables linéaires et non linéaires

Que nous montrent ces graphiques sur la relation entre la bioluminescence, la profondeur de l'échantillon et les saisons?

```r
par(mfrow = c(1,2))
plot(basic_model, all.terms = TRUE)
```

???

La bioluminescence varie de manière non-linéaire sur le gradient `SampleDepth`, avec des niveaux de bioluminescence les plus élevés à la surface, suivis d'un second maximum plus petit, juste au-dessus d'une profondeur de 1500, avec des niveaux décroissants à des profondeurs plus basses.

Il y a également une différence prononcée dans la bioluminescence entre les saisons, avec des niveaux élevés pendant la saison 2, par rapport à la saison 1.

---
# Degrés de liberté effectifs (EDF)

```r
basic_summary$s.table
#                     edf   Ref.df        F p-value
# s(SampleDepth) 8.706426 8.975172 184.3583       0
```

Les `edf` indiqués dans le `s.table` sont les **degrés effectifs de liberté** (EDF, en anglais "_effective degrees of freedom_") du terme lisse `s(SampleDepth)`.

Plus le nombre de degrés de liberté est élevé, plus la courbe est complexe et ondulée.

- Une valeur proche de 1 se rapproche d'un terme linéaire.

- Une valeur elevée signifie que la courbe est plus ondulé, ou en d'autres termes, plus non-linéaire.

> Dans notre modèle de base, les **EDF** du terme non-linéaire `s(SampleDepth)` sont ~9, ce qui suggère une courbe fortement non-linéaire.

---
# Degrés de liberté effectifs (EDF)

Les degrés de liberté effectifs nous donnent beaucoup d'informations sur la relation entre les prédicteurs du modèle et les variables de réponse.

> Vous reconnaissez peut-être le terme "degrés de liberté" suite à des ateliers précédents sur les modèles linéaires, **mais attention**!

> Les degrés de liberté effectifs d'un GAM sont estimés différemment des degrés de liberté d'une régression linéaire, et sont interprétés différemment.

Dans la régression linéaire, les degrés de liberté du *modèle* sont équivalents au nombre de paramètres libres non redondants, `$p$`, dans le modèle (et les degrés de liberté *résiduels* sont égaux à `$n-p$`; où `$n$` est le nombre d'observations).

Parce que le nombre de paramètres libres des fonctions de lissage est souvent difficile à définir, les **EDF** sont liés à `$\lambda$`, où l'effet de la pénalité est de réduire les degrés de liberté.

???

Revenons sur le concept de degrés de liberté effectifs (EDF). (Rappel: Défi 1)

---

# Degrés de liberté effectifs (EDF) et `$k$`

La limite supérieure d'**EDF** est déterminée par les dimensions de base `$k$` de la fonction lisse (les **EDF** ne peuvent pas dépasser `$k-1$`)

En pratique, le choix exact de `$k$` est arbitraire, mais il devrait être **suffisamment grand** pour permettre une fonction lisse suffisamment complexe.

Nous discuterons du choix de `$k$` dans les sections qui suivent.

---

# GAM à plusieurs variables linéaires et lisses

Nous pouvons ajouter un second terme, `RelativeDepth`, mais en spécifiant une relation linéaire avec `Sources`:

```r
two_term_model <- gam(Sources ~ Season + s(SampleDepth) + RelativeDepth,
 data = isit, method = "REML")
two_term_summary <- summary(two_term_model)
```

Informations sur les effets paramétriques (termes linéaires):

```r
two_term_summary$p.table
#                   Estimate   Std. Error   t value     Pr(>|t|)
# (Intercept)    9.808305503 0.6478741951 15.139213 1.446613e-45
# Season2        6.041930627 0.4767977508 12.671894 1.380010e-33
# RelativeDepth -0.001401908 0.0002968443 -4.722705 2.761048e-06
```

Informations sur les effets additifs (termes non linéaires):

```r
two_term_summary$s.table
#                     edf  Ref.df        F p-value
# s(SampleDepth) 8.699146 8.97396 132.4801       0
```

???

L'estimation du coefficient de régression pour ce nouveau terme linéaire, `RelativeDepth`, sera présenté dans le tableau `p.table`. Rappelez-vous, le tableau `p.table` montre des informations sur les effets paramétriques (termes linéaires).

Dans `s.table`, nous trouverons encore une fois le terme non-linéaire, `s(SampleDepth)`, et son paramètre de courbure (`edf`). Rappelez-vous, le tableau `s.table` montre des informations sur les effets additifs (termes non-linéaires).

---

# GAM à plusieurs variables linéaires et lisses

```r
par(mfrow = c(2,2))
plot(two_term_model, all.terms = TRUE)
```

---

# GAM à plusieurs variables non linéaires

Si nous voulons vérifier que la relation entre `Sources` et `RelativeDepth` est non-linéaire, on peut modéliser `RelativeDepth` avec une fonction non-linéaire.

```r
two_smooth_model <- gam(Sources ~ Season + s(SampleDepth) + s(RelativeDepth),
 data = isit, method = "REML")
two_smooth_summary <- summary(two_smooth_model)
```

Informations sur les effets paramétriques (termes linéaires) :

```r
two_smooth_summary$p.table
#             Estimate Std. Error  t value     Pr(>|t|)
# (Intercept) 7.937755  0.3452945 22.98836 1.888513e-89
# Season2     4.963951  0.4782280 10.37988 1.029016e-23
```

Informations sur les effets additifs (termes non linéaires) :

```r
two_smooth_summary$s.table
#                       edf   Ref.df         F p-value
# s(SampleDepth)   8.752103 8.973459 150.37263       0
# s(RelativeDepth) 8.044197 8.749580  19.97476       0
```

???

---
# GAM à plusieurs variables non linéaires

Nous pouvons aussi vérifier que la relation entre `Sources` et `RelativeDepth` est non-linéaire.

```r
par(mfrow = c(2,2))
plot(two_smooth_model, all.terms = TRUE)
```

???

Pensez-vous que la performance de notre modèle est amélioré par l'ajout de ce nouveau terme non-linéaire, pour mieux représenter la relation entre la bioluminescence et la profondeur relative?

---
# GAM à plusieurs variables non linéaires

Comme précédemment, nous pouvons comparer nos modèles avec AIC pour tester si le terme non-linéaire améliore la performance de notre modèle:

```r
AIC(basic_model, 
    two_term_model, 
    two_smooth_model)
#                        df      AIC
# basic_model      11.83374 5208.713
# two_term_model   12.82932 5188.780
# two_smooth_model 20.46960 5056.841
```

On peut voir que `two_smooth_model` a la plus petite valeur de AIC.

]

```r
two_smooth_model
# 
# Family: gaussian 
# Link function: identity 
# 
# Formula:
# Sources ~ Season + s(SampleDepth) + s(RelativeDepth)
# 
# Estimated degrees of freedom:
# 8.75 8.04  total = 18.8 
# 
# REML score: 2550.112
```

Le modèle le mieux ajusté comprend donc deux fonctions non-linéaires pour `SampleDepth` et `RelativeDepth`, et un terme linéaire pour `Season`.

]

---
exclude:true
# GAM à plusieurs variables

Nous pouvons ajouter un second terme, `x2`, mais spécifier une relation linéaire avec `y`

```r
two_term_model <- gam(y ~ x0 + s(x1) + x2, data = gam_data)
two_term_summary <- summary(two_term_model)
two_term_summary$p.table

two_term_summary$s.table
```

---
exclude:true
# GAM à plusieurs variables

Nous pouvons ajouter un second terme, `x2`, mais spécifier une relation linéaire avec `y`

```r
plot(two_term_model, all.terms = TRUE)
```

---
exclude:true
# GAM à plusieurs variables

Nous pouvons aussi explorer si la relation entre `y` et `x2` est non-linéaire

```r
two_smooth_model <- gam(y ~ x0 + s(x1) + s(x2), data = gam_data)
two_smooth_summary <- summary(two_smooth_model)
two_smooth_summary$p.table

two_smooth_summary$s.table
```

---
exclude:true
# GAM à plusieurs variables

Nous pouvons aussi explorer si la relation entre `y` et `x2` est non-linéaire

```r
plot(two_smooth_model, page = 1, all.terms = TRUE)
```

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

Pour notre deuxième défi, nous allons développer notre modèle en ajoutant des variables qui, selon nous, pourraient être des prédicteurs écologiquement significatifs pour expliquer la bioluminescence.

1. Créez deux nouveaux modèles: Ajoutez la variable `Latitude` à `two_smooth_model`, premièrement comme paramètre linéaire, et ensuite comme fonction non-linéaire.

2. Est-ce que `Latitude` est un terme important à inclure dans le modèle? La `Latitude` a-t-elle un effet linéaire ou non-linéaire?

> Utilisez des graphiques, les tables des coefficients et la fonction `AIC()` pour répondre à ces questions.

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

1. Créez deux nouveaux modèles avec la variable `x3` comme paramètre linéaire et non linéaire.
2. Utilisez des graphiques, les tables des coefficients et la fonction `AIC()` afin de déterminer s'il est nécessaire d'inclure `x3` dans le modèle.

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

**1.** Créez deux nouveaux modèles: Ajoutez la variable `Latitude` à `two_smooth_model`, premièrement comme paramètre linéaire, et ensuite comme fonction non-linéaire.

```r
# Ajouter Latitude comme terme linéaire
three_term_model <- gam(Sources ~ 
 Season + s(SampleDepth) + s(RelativeDepth) + 
* Latitude,
 data = isit, method = "REML")

three_term_summary <- summary(three_term_model)
```

```r
# Ajouter Latitude comme terme non-linéaire
three_smooth_model <- gam(Sources ~ 
 Season + s(SampleDepth) + s(RelativeDepth) + 
* s(Latitude),
 data = isit, method = "REML")

three_smooth_summary <- summary(three_smooth_model)
```

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

__2.__ Est-ce que `Latitude` est un terme important à inclure dans le modèle?

Commençons par visualiser les 4 effets qui sont maintenant inclus dans chaque modèle.

```r
par(mfrow = c(2,2))

plot(three_smooth_model, 
     all.terms = TRUE)
```
]

]

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

```r
three_term_model <- gam(y ~ x0 + s(x1) + s(x2) + x3, data = gam_data)
three_smooth_model <- gam(y~x0 + s(x1) + s(x2) + s(x3), data = gam_data)
three_smooth_summary <- summary(three_smooth_model)
```

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

```r
par(mfrow = c(2,2))
plot(three_smooth_model, all.terms = TRUE)
```

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

Nous devrions également examiner nos tableaux de coefficients. Qu'est-ce que les EDF nous disent à propos de _l'ondulation_, ou la non-linéarité, des effets de nos prédicteurs?

```r
three_smooth_summary$s.table
#                       edf   Ref.df         F p-value
# s(SampleDepth)   8.766891 8.975682 68.950905       0
# s(RelativeDepth) 8.007411 8.730625 17.639321       0
# s(Latitude)      7.431116 8.296838  8.954349       0
```

Les EDF sont tous élevés pour nos variables, y compris `Latitude`.

Cela nous indique que `Latitude` est assez _ondulée_, et qu'elle ne devrait probablement pas être incluse comme terme linéaire.

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

Avant de décider quel modèle est le "meilleur", nous devrions tester si l'effet `Latitude` est plus approprié comme terme linéaire ou lisse, en utilisant `AIC()`:

```r
AIC(three_smooth_model, three_term_model)
#                          df      AIC
# three_smooth_model 28.20032 4990.546
# three_term_model   21.47683 5051.415
```

Notre modèle incluant la `Latitude` comme terme _non-linéaire_ a un score AIC inférieur, ce qui signifie qu'il est plus performant que notre modèle incluant la `Latitude` comme terme _linéaire_.

Mais, est-ce que l'ajout de `Latitude` comme prédicteur non-linéaire améliore réellement notre "meilleur" modèle (`two_smooth_model`)?

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

```r
AIC(two_smooth_model, three_smooth_model)
#                          df      AIC
# two_smooth_model   20.46960 5056.841
# three_smooth_model 28.20032 4990.546
```

Notre `three_smooth_model` a un score AIC inférieur à notre meilleur modèle précédent (`two_smooth_model`), qui n'incluait pas `Latitude`.

Ceci implique que `Latitude` est en effet un prédicteur informatif non-linéaire de la bioluminescence.

---
class: inverse, center, middle

## 4. Interactions

---
# GAM avec des termes d'interaction

Il y a deux façons de modéliser une interaction entre deux variables :

- pour __deux variables non-linéaires__ : `s(x1, x2)`

- pour __une variable non-linéaire et une variable linéaire__ (quantitative ou qualitative) : utiliser l'argument `by`, `s(x1, by = x2)`
- Quand `x2` est qualitative, vous avez un terme non linéaire qui varie entre les différents niveaux de `x2`
- Quand `x2` est quantitative, l'effet linéaire de `x2` varie avec `x1`
- Quand `x2` est qualitative, le facteur doit être ajouté comme effet principal dans le modèle

---
# Interaction: variables non-linéaire et qualitatif

Nous allons examiner l'effet de l'interaction en utilisant notre variable qualitative `Season` et examiner si la non-linéarité de `s(SampleDepth)` varie selon les différents niveaux de `Season`.

```r
factor_interact <- gam(Sources ~ Season +
 s(SampleDepth,by=Season) +
 s(RelativeDepth),
 data = isit, method = "REML")
```

```r
summary(factor_interact)$s.table
#                             edf   Ref.df          F p-value
# s(SampleDepth):Season1 6.839386 7.552045  95.119422       0
# s(SampleDepth):Season2 8.744574 8.966290 154.474325       0
# s(RelativeDepth)       6.987223 8.055898   6.821074       0
```

---
exclude:true
# GAM avec des termes d'interaction

Nous allons examiner l'effet de l'interaction en utilisant notre variable qualitative `x0` et examiner si la non-linéarité de `s(x2)` varie selon les différents niveaux de `x0`.

```r
factor_interact <- gam(y ~ x0 + s(x1) + s(x2, by = x0), data = gam_data)

summary(factor_interact)$s.table
```

---
# Interaction: variables non-linéaire et qualitatif

```r
par(mfrow = c(2),2))
plot(factor_interact)
```

???

Les deux graphiques montrent l'effet d'interaction entre notre variable lisse `SampleDepth` et chaque niveau de notre variable factorielle, `Season`.

Une bonne question pour les participants: Voyez-vous une différence entre les deux courbes?

Les graphiques montrent quelques différences entre la forme des termes lisses entre les deux niveaux de `Season`. La différence la plus notable est le pic dans le deuxième panneau, qui nous indique qu'il y a un effet de `SampleDepth` entre 1000 et 2000 qui est important dans la saison 2, mais qui ne se produit pas dans la saison 1.

Ceci suggère que l'effet d'interaction pourrait être important à inclure dans notre modèle.

---
# Interaction: variables non-linéaire et qualitatif

Nous pouvons également représenter l'effet d'interaction en 3D sur un seul graphique, en utilisant `vis.gam()`.

```r
vis.gam(factor_interact, theta = 120, 
        n.grid = 50, lwd = .4)
```

]

```r
vis.gam(factor_interact, theta = 60, 
        n.grid = 50, lwd = .4)
```

]

???

Nous pouvons modifier la rotation de ce graphique en utilisant l'argument `theta`.

Lorsqu'on visualise les termes d'interaction, on observe des différences dans la forme du lissage de `SampleDepth` entre les deux niveaux de `Season`. Ceci suggère que l'effet d'interaction pourrait être important à inclure dans notre modèle.

---
# Interaction: variables non-linéaire et qualitatif

Ces graphiques suggèrent que l'effet d'interaction pourrait être important à inclure dans notre modèle.

On peut faire une comparaison de modèles en utilisant l'AIC pour déterminer si le terme d'interaction améliore la performance de notre modèle.

```r
AIC(two_smooth_model, factor_interact)
#                        df      AIC
# two_smooth_model 20.46960 5056.841
# factor_interact  26.99693 4878.631
```

L'AIC de notre modèle avec une interaction factorielle entre la variable non-linéair `SampleDepth` et le `Season` a un score AIC plus bas.

> Ensemble, l'AIC et les graphiques confirment que l'inclusion de cette interaction améliore la performance de notre modèle.

---
# Interaction entre variables non linéaires

Finalement, nous regardons les interactions entre deux termes non linéaires, `SampleDepth` et `RelativeDepth`.

```r
smooth_interact <- gam(Sources ~ Season + s(SampleDepth, RelativeDepth),
 data = isit, method = "REML")
summary(smooth_interact)$s.table
# edf Ref.df F p-value
# s(SampleDepth,RelativeDepth) 27.12521 28.77 93.91722 0
```

---

# Interaction entre variables non linéaires

```r
vis.gam(smooth_interact,  view = c("SampleDepth", "RelativeDepth"), 
        theta = 50, n.grid = 50, lwd = .4)
```

???

Les graphiques illustrent une interaction non linéaire, où `Sources` est plus faible à des valeurs élevées de `SampleDepth` et `RelativeDepth`, mais augmente avec `RelativeDepth` alors que `SampleDepth` est faible.

Rappelez-vous, ce graphique peut être réorienté en changeant la valeur de l'argument `theta`.
On peut changer la couleur du graphique 3D en utilisant l'argument `color`. Essayez de spécifier `color = "cm"` dans `vis.gam()` ci-dessus, et consultez `?vis.gam` pour plus d'options de couleurs.

---
# Interaction entre variables non linéaires

Ainsi, il semble y avoir un effet d'interaction entre ces termes non linéaires.

Est-ce que l'inclusion de l'interaction entre `s(SampleDepth)` et `s(RelativeDepth)` améliore notre modèle `two_smooth_model`?

```r
AIC(two_smooth_model, smooth_interact)
#                        df      AIC
# two_smooth_model 20.46960 5056.841
# smooth_interact  30.33625 4943.890
```

Le modèle avec l'intéraction entre `s(SampleDepth)` et `s(RelativeDepth)` a une plus petite valeur d'AIC.

> L'inclusion de cette interaction améliore la performance de notre modèle, et donc notre capacité à comprendre les déterminants de la bioluminescence.

---
class: inverse, center, middle

# 5. Généralisation du modèle additif

---

# Généralisation du modèle additif

Le modèle additif de base peut être étendu de plusieurs façons :

1. Utiliser d'autres __distributions__ pour la variable de réponse avec l'argument `family` (comme dans un GLM),
2. Utiliser de différents types de __fonctions de base__,
3. Utilisation de différents types d'__effets aléatoires__ pour ajuster des modèles à effets mixtes.

Nous allons maintenant examiner ces aspects.

---

# Modèle additif généralisé

Jusqu'à présent, nous avons utilisé des modèles additifs simples (gaussiens), l'équivalent non linéaire d'un modèle linéaire.

Mais que pouvons-nous faire si :
- les observations de la variable de réponse ne **suivent pas une distribution Normale**?
- la **variance n'est pas constante**? (hétéroscédasticité)

Tout comme les modèles linéaires généralisés (GLM), nous pouvons formuler des modèles additifs **généralisés** pour répondre à ces problèmes.

---

# Modèle additif généralisé

Revenons au modèle d'interaction pour les données de bioluminescence :

```r
smooth_interact <- gam(Sources ~ Season + s(SampleDepth, RelativeDepth),
 data = isit, method = "REML")

summary(smooth_interact)$p.table
#             Estimate Std. Error   t value     Pr(>|t|)
# (Intercept) 8.077356  0.4235432 19.070912 1.475953e-66
# Season2     4.720806  0.6559436  7.196969 1.480113e-12

summary(smooth_interact)$s.table
#                                   edf Ref.df        F p-value
# s(SampleDepth,RelativeDepth) 27.12521  28.77 93.91722       0
```

---

# Validation d'un GAM

Comme pour un GLM, il est essentiel de vérifier si nous avons correctement spécifié le modèle, en particulier la *distribution* de la réponse.

Il faut vérifier:

1. Le choix des dimensions de base `k`.
2. Les tracés des résidus (comme pour un GLM).

Fonctions incluses dans `mgcv` :

- `k.check()` effectue une vérification des dimensions de base.
- `gam.check()` produit des tracés de résidus (et fournit également la sortie de `k.check()`.

---
# Validation GAM: choisir `$k$` dimensions de base

Vous vous rappelez du paramètre de lissage `$\lambda$` qui restreint les _ondulations_ de nos fonctions de lissage?

Cette _ondulation_ est également contrôlé par la dimension de base `$k$`, qui définit le nombre de fonctions de base utilisées pour créer une fonction lisse.

Plus le nombre de fonctions de base utilisées pour construire une fonction lisse est élevé, plus la fonction lisse est "ondulée":

???

Dans les section précédents, nous avons discuté du rôle du paramètre de lissage `$\lambda$` pour contrôler les _ondulations_ de nos fonctions de lissage. Cette _ondulation_ est également contrôlé par la dimension de base `$k$`, qui définit le nombre de fonctions de base utilisées pour créer une fonction lisse.

Chaque fonction lisse dans un GAM est essentiellement la somme pondérée de nombreuses fonctions plus petites, appelées fonctions de base. Plus le nombre de fonctions de base utilisées pour construire une fonction lisse est élevé, plus la fonction lisse est "ondulée". Comme vous pouvez le voir ci-dessous, une fonction lisse avec une petite dimension de base de `$k$` sera moins ondulée qu'une fonction lisse avec une grande dimension de base de `$k$`.

---
# Validation GAM: choisir `$k$` dimensions de base

La clé pour obtenir un bon ajustement du modèle consiste à __équilibrer__ le compromis entre deux éléments:

+ Le paramètre de lissage `$\lambda$`, qui _pénalise les ondulations_ ;
+ La dimension de base `$k$`, qui permet plus de _flexibilité_ au modèle (plus d'ondulations) en fonction de nos données.

Avons-nous optimisé le compromis entre le lissage ( `$\lambda$` ) et la flexibilité ( `$k$` ) dans notre modèle?

En d'autres mots, est-ce que le `k` est assez grand?

```r
k.check(smooth_interact)
#                              k'      edf   k-index p-value
# s(SampleDepth,RelativeDepth) 29 27.12521 0.9448883  0.0575
```

Les **EDF se rapprochent beaucoup de** `k`, donc la _flexibilité_ du modèle est trop restreinte par le `k` utilisé par défaut.

En d'autres mots, .alert[nous n'avons pas un compromis équilibré entre le lissage et l'ondulation du modèle].

???

Ici, on se demande essentiellement: Est-ce que notre module est assez flexible?

On n'a pas encore spécifié de valeur `$k$` dans notre modèle, mais `gam()` définit un `$k$` par défaut en fonction du nombre de variables sur lequel la fonction lisse est construite.

Les **EDF se rapprochent beaucoup de** `k`. Ceci signifie que la _flexibilité_ du modèle est restreint par le `k` utilisé par défaut, et que notre modèle pourrait mieux s'ajuster aux données si on permettait plus d'ondulations. En d'autres mots, nousn n'avons pas un compromis équilibré entre le lissage et l'ondulation du modèle.

---
### Validation GAM: choisir `$k$` dimensions de base

Recommençons le modèle avec un `k` plus élevé:

```r
smooth_interact_k60 <- gam(Sources ~ Season + s(SampleDepth, RelativeDepth, k = 60),
 data = isit, method = "REML")
```

```r
summary(smooth_interact_k60)$p.table
#             Estimate Std. Error   t value     Pr(>|t|)
# (Intercept) 8.129512  0.4177659 19.459491 1.911267e-68
# Season2     4.629964  0.6512205  7.109672 2.741865e-12
summary(smooth_interact_k60)$s.table
#                                   edf   Ref.df        F p-value
# s(SampleDepth,RelativeDepth) 46.03868 54.21371 55.19817       0
```

Est-ce que `k` est assez grand maintenant?

```r
k.check(smooth_interact_k60)
#                              k'      edf  k-index p-value
# s(SampleDepth,RelativeDepth) 59 46.03868 1.048626   0.895
```

Les EDF sont beaucoup plus petits que `k`. On peut remplacer notre modèle avec cette version plus flexible:

```r
smooth_interact <- smooth_interact_k60
```

???

Les EDF sont beaucoup plus petits que `k`, donc notre modèle s'adjuste mieux aux données avec plus d'ondulations. On peut donc remplacer notre modèle avec cette version plus flexible.

---
# Validation GAM: Choisir une autre distribution

Comme pour tout modèle normal, nous devons vérifier certaines conditions avant de continuer.

```r
par(mfrow = c(2,2))
gam.check(smooth_interact)
```

]

]

???

Rappel: Nous pouvons évaluer la distribution des résidus du modèle pour vérifier ces conditions, tout comme nous le ferions pour un GLM (voir [Atelier 6](https://r.qcbs.ca/fr/workshops/r-workshop-06/)).

---
# Validation GAM: Choisir une autre distribution

--

???

- Ces tracés sont un peu différents de ceux produits par `plot` pour un modèle linéaire (par exemple, pas de tracé de levier).

- Les participants devraient déjà être familiarisés avec les graphiques de résidus (ils sont expliqués plus en détail dans les [Atelier 4](https://r.qcbs.ca/fr/workshops/r-workshop-04/) et [Atelier 6](https://r.qcbs.ca/fr/workshops/r-workshop-06/).

La visualisation des résidus met quelques problèmes en évidence:
- Figure 2: La variance des erreurs n'est pas constante (hétéroscédasticité).
- Figures 1 et 4: Quelques observations extrêmes.

---
# Validation GAM: Choisir une autre distribution

Pour notre modèle d'interaction, nous avons besoin d'une distribution de probabilité qui permet à la **variance d'augmenter avec la moyenne**.

Une famille de distributions qui possède cette propriété et qui fonctionne bien dans un GAM est la famille **Tweedie**.

Une fonction de liaison commune pour les distributions *Tweedie* est le `$log$`.

Comme dans un GLM, nous pouvons utiliser l'argument `family = ` dans `gam()` pour ajuster des modèles avec d'autres distributions (y compris des distributions telles que `binomial`, `poisson`, `gamma` etc).

Pour en savoir plus sur les familles disponibles dans `mgcv` :

```r
?family.mgcv
```

---

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

1. Ajuster un nouveau modèle `smooth_interact_tw` avec la même formule que le modèle `smooth_interact` mais avec une distribution de la famille *Tweedie* (au lieu de la distribution normale) et `log` comme fonction de liaison. Pour ce faire, on peut utiliser `family = tw(link = "log")` dans `gam()`.
2. Vérifier le choix de `k` et les tracés de résidus pour le nouveau modèle.
3. Comparer `smooth_interact_tw` avec `smooth_interact`. Lequel est meilleur ?

```r
# Voici comment nous écririons le modèle pour spécifier la distribution normale :
smooth_interact <- gam(Sources ~ Season + s(SampleDepth, RelativeDepth, k = 60),
 family = gaussian(link = "identity"),
 data = isit, method = "REML")
```

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

__1.__ Premièrement, construisons un nouveau modèle avec une distribution *Tweedie* un lien `log`.

```r
smooth_interact_tw <- gam(Sources ~ Season + s(SampleDepth, RelativeDepth, k = 60),
 family = tw(link = "log"),
 data = isit, method = "REML")
summary(smooth_interact_tw)$p.table
# Estimate Std. Error t value Pr(>|t|)
# (Intercept) 1.3126641 0.03400390 38.60333 8.446478e-180
# Season2 0.5350529 0.04837342 11.06089 1.961733e-26
summary(smooth_interact_tw)$s.table
# edf Ref.df F p-value
# s(SampleDepth,RelativeDepth) 43.23949 51.57139 116.9236 0
```

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

__2.__ Vérifier le choix de `k` et la visualisation des résidus pour le nouveau modèle.

Ensuite, on devrait vérifier le choix des dimensions de bases `k`:

```r
k.check(smooth_interact_tw)
#                              k'      edf  k-index p-value
# s(SampleDepth,RelativeDepth) 59 43.23949 1.015062  0.7825
```

On a un équilibre entre `$k$` et les EDF. Super!

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

Ensuite la visualisation des résidus, pour valider la distribution Tweedie:

```r
par(mfrow=c(2,2))
gam.check(smooth_interact_tw)
```

???
Les résidus semblent mieux distribués, mais il est évident que le modèle manque encore quelque chose. Il pourrait s'agir d'un effet spatial (longitude et latitude), ou d'un effet aléatoire (par exemple basé sur `Station`).

---

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

__3.__ Comparer `smooth_interact_tw` avec `smooth_interact`. Lequel est meilleur?

```r
AIC(smooth_interact, smooth_interact_tw)
#                          df      AIC
# smooth_interact    49.47221 4900.567
# smooth_interact_tw 47.86913 3498.490
```

.comment[L'AIC nous permet de comparer des modèles qui sont basés sur des distributions différentes!]

Utiliser une distribution *Tweedie* au lieu d'une distribution *Normale* améliore beaucoup notre modèle!

???

Le score AIC pour `smooth_interact_tw` est beaucoup plus petit que celui de `smooth_interact`. Utiliser une distribution *Tweedie* au lieu d'une distribution *Normale* améliore beaucoup notre modèle!

---
class: inverse, center, middle

## 6. Changer la fonction de base

---
# Autres fonctions lisses

Pour modéliser une surface lisse ou non-linéaire, nous pouvons construire des fonctions lisses de différentes manières:

`s()` ![:faic](arrow-right) pour modéliser un terme lisse 1-dimensionnelle, ou pour modéliser une intéraction entre des variables mesurées sur la *même échelle*

`te()` ![:faic](arrow-right)  pour modéliser une surface d'interaction 2- ou n-dimensionnelle entre des variables qui *ne sont pas sur la même échelle*. Comprend les effets principaux.

`ti()` ![:faic](arrow-right) pour modéliser une surface d'interaction 2- ou n-dimensionnelle *qui ne comprend pas les effets principaux*.

---
# Paramètres des fonctions lisses

Les fonctions lisses ont beaucoup de paramètres qui pourraient changer leur comportement. Les paramètres les plus souvent utilisés sont les suivants :

`k` ![:faic](arrow-right) dimensions de base
- détermine la limite supérieure du nombre de fonctions de base utilisées pour construire la courbe.
- contraint l'ondulation d'une fonction lisse.
- k devrait être < au nombre de données uniques.
- la complexité (ou la non-linéarité) d'une fonction lisse dans un modèle ajusté est reflétée par ses degrés de liberté effectifs (**EDF**)

---
# Paramètres des fonctions lisses

Les fonctions lisses ont beaucoup de paramètres qui pourraient changer leur comportement. Les paramètres les plus souvent utilisés sont les suivants :

`bs` ![:faic](arrow-right) spécifie la fonction de base sous-jacente.
- pour `s()` on utilise `tp` (*thin plate regression spline*) et pour `te()` et `ti()` on utilise la base `cr` (*cubic regression spline*).

`d` ![:faic](arrow-right) spécifie quelles variables d'une intéraction se trouvent sur la même échelle lorsqu'on utilise `te()` and `ti()`.
- Par exemple, `te(Temps, largeur, hauteur, d=c(1,2))`, indique que la `largeur` et la `hauteur` sont sur la même échelle, mais que `temps` ne l'est pas.

---
# Exemple: Données cycliques

Lorsque l'on modélise des données cycliques, on souhaite généralement que le prédicteur soit identique aux deux bouts des phases. Pour y parvenir, nous devons modifier la fonction de base.

Utilisons une série chronologique de données climatiques, divisées en mesures mensuelles, afin de déterminer s'il y a une tendance de température annuelle. Nous utiliserons la série temporelle de la température de Nottingham dans `R`:

```r
data(nottem)
```

```r
# nombre d'années de données (20 ans)
n_years <- length(nottem)/12

# codage qualitatif pour les 12 mois de l'année,
# pour chaque année échantillonnée (série de 1 à 12, répétée 20 fois)
nottem_month <- rep(1:12, times = n_years)

# une variable où l'année correspondant à chaque mois dans nottem_month
nottem_year <- rep(1920:(1920 + n_years - 1), each = 12)
```

```r
# Visualiser la série temporelle
qplot(x = nottem_month, y = nottem, colour = factor(nottem_year), geom = "line") +
  theme_bw()
```

???

Voir `?nottem` pour plus de détails sur le jeu de données.

---
# Exemple: Données cycliques

---
# Exemple: Données cycliques

Nous pouvons modéliser le changement cyclique de température à travers les mois et la tendance non-linéaire à travers les années, en utilisant une spline cubique, ou `cc` pour modéliser les effets de mois ainsi qu'un terme non-linéaire pour la variable année.

```r
year_gam <- gam(nottem ~ s(nottem_year) + s(nottem_month, bs = "cc"), 
 method = "REML")
```

```r
summary(year_gam)$s.table
#                      edf   Ref.df          F    p-value
# s(nottem_year)  1.621375 2.011475   2.788093 0.06141004
# s(nottem_month) 6.855132 8.000000 393.119285 0.00000000
```

---
# Exemple: Données cycliques

```r
plot(year_gam, page = 1, scale = 0)
```

Il y a une augmentation de 1 à 1.5 degrés au cours de la série, mais au cours d'une année, il y a une variation d'environ 20 degrés. Les données réelles varient autour de ces valeurs prédites et ceci représente donc la variance inexpliquée.

???

Ici, nous pouvons voir l'un des avantages très intéressants de l'utilisation des GAMs. Nous pouvons soit tracer la surface réponse (valeurs prédites) ou les termes (contribution de chaque covariable) tel qu'indiqué ci-haut. Vous pouvez imaginer ce dernier en tant qu'une illustration de la variation des coefficients de régression et comment leur contribution (ou taille de leur effet) varie au fil du temps. Dans le premier graphique, nous voyons que les contributions positives de la température sont survenues après 1930.

---
class: inverse, center, middle

# 7. Introduction rapide aux GAMMs

---

# La non-indépendance des données

Lorsque les observations ne sont pas indépendantes, les GAMs peuvent
être utilisés soit pour incorporer:

- une structure de corrélation pour modéliser les résidus autocorrélés, comme:
 - le modèle autorégressif (AR, en anglais _autoregressive_); 
 - le modèle avec moyenne mobile (MA, en anglais _moving average_); ou,
 - une combinaison des deux modèles (ARMA, en anglais, _autoregressive moving average_).
- des effets aléatoires qui modélisent l'indépendance entre les observations d'un même site.

---
# L'autocorrelation des résidus

**L'autocorrélation des résidus** désigne le degré de corrélation entre les résidus (les différences entre les valeurs réelles et les valeurs prédites) d'un modèle de série temporelle.

En d'autres termes, s'il y a autocorrélation des résidus dans un modèle de série temporelle, cela signifie qu'il existe un modèle ou une relation entre les résidus à un moment donné et les résidus à d'autres moments.

L'autocorrélation des résidus est généralement mesurée à l'aide des graphiques **ACF (fonction d'autocorrélation)** et **pACF (fonction d'autocorrélation partielle)**, qui montrent la corrélation entre les résidus à différents retards.

Si les graphiques **ACF** ou **pACF** montrent des corrélations significatives à des retards non nuls, il y a des preuves d'autocorrélation des résidus et le modèle peut devoir être modifié ou amélioré pour mieux capturer les modèles sous-jacents dans les données.

Voyons comment cela fonctionne avec notre modèle `year_gam`!

---
# L'autocorrelation des résidus

Pour commencer, nous allons jeter un coup d'œil à un modèle avec de l'autocorrélation temporelle dans les résidus.

Revenons au modèle de la température de Nottingham pour vérifier si les résidus sont corrélés en faisant appel à la fonction (partielle) d'autocorrélation.

```r
par(mfrow = c(1,2))

acf(resid(year_gam), lag.max = 36, main = "ACF")

pacf(resid(year_gam), lag.max = 36, main = "pACF")
```

Souvenez vous que: `acf` signifie _autocorrelation function_, et `pacf` signifie _partial autocorrelation function_.

---
# L'autocorrelation des résidus

.comment[__ACF__ évalue la corrélation croisée et __pACF__ la corrélation partielle d'une série temporelle entre points à différents décalages (donc, la similarité entre observations progressivement décalés).]

Les graphiques ACF et pACF sont donc utilisés pour identifier le temps nécessaire avant que les observations __ne sont plus autocorrélées__.

Un modèle AR de faible ordre semble nécessaire (donc, avec un ou deux intervalles de temps décalés).

???

La __fonction d'autocorrelaton__ (ACF; première figure) évalue la corrélation croisée d'une série temporelle entre points à différents décalages (donc, la similarité entre observations progressivement décalés).

En contraste, la __fonction partielle d'autocorrelation__ (PACF: deuxième figure) évalue la corrélation croisée d'une série temporelle entre points à différents décalages, après avoir contrôlé les valeurs de la série temporelle à tous les décalages plus courts.

Dans le graphique ACF, la région ombrée en bleu représente l'intervalle de confiance et les lignes rouges en pointillé représentent les limites de la signification statistique.

---
# L'autocorrelation des résidus

Ajoutons des structures d'autocorrelation au modèle de température de Nottingham:

- **AR(1)** : corrélation avec un intervalle de temps décalé, ou
- **AR(2)** : corrélation à deux intervalles de temps décalés.

```r
df <- data.frame(nottem, nottem_year, nottem_month)

year_gam <- gamm(nottem ~ s(nottem_year) + s(nottem_month, bs = "cc"), data = df)

year_gam_AR1 <- gamm(nottem ~ s(nottem_year) + s(nottem_month, bs = "cc"),
 correlation = corARMA(form = ~ 1|nottem_year, p = 1),
 data = df)

year_gam_AR2 <- gamm(nottem ~ s(nottem_year) + s(nottem_month, bs = "cc"),
 correlation = corARMA(form = ~ 1|nottem_year, p = 2),
 data = df)
```

---

# L'autocorrelation des résidus

Quel modèle est mieux ajusté?

```r
AIC(year_gam$lme, year_gam_AR1$lme, year_gam_AR2$lme)
#                  df      AIC
# year_gam$lme      5 1109.908
# year_gam_AR1$lme  6 1101.218
# year_gam_AR2$lme  7 1101.598
```

Le modèle avec la structure AR(1) donne un meilleur ajustement comparé au premier modèle (`year_gam`), il y a très peu d'amélioration en passant au `AR(2)`.

Il est donc préférable d'inclure uniquement la structure `AR(1)` dans notre modèle.

---
# Effets aléatoires

Comme nous l'avons vu dans la section précédente, `bs` spécifie la fonction de base sous-jacente. Pour les facteurs aléatoires (origine et pente linéaire), nous utilisons `bs = "re"` et pour les pentes aléatoires non linéaires, nous utilisons `bs = "fs"`.

--

**Trois types d'effets aléatoires différents** sont possibles lors de l'utilisation des GAMMs (où `fac` ![:faic](arrow-right) variable qualitative utilisée pour l'effet aléatoire; `x0` ![:faic](arrow-right) effet quantitatif fixe) :

- **interceptes aléatoires** ajustent la hauteur des termes du modèle avec une valeur constante de pente : `s(fac, bs = "re")`

- **pentes aléatoires** ajustent la pente d'une variable explicative numérique : `s(fac, x0, bs = "re")`

- **surfaces lisses aléatoires** ajustent la tendance d'une prédiction numérique de façon non linéaire: `s(x0, fac, bs = "fs", m = 1)`, où l'argument `$m = 1$` met une plus grande pénalité au lissage qui s'éloigne de 0, ce qui entraîne un retrait vers la moyenne.

---
# GAMM avec un intercepte aléatoire

Nous allons utiliser `gamSim()` pour simuler un ensemble de données, cette fois-ci avec un effet aléatoire. Ensuite, nous construirons un modèle avec un intercepte aléatoire en utilisant `fac` comme facteur aléatoire.

```r
# Simuler des données
gam_data2 <- gamSim(eg = 6)

head(gam_data2)
```

```
# 4 term additive + random effectGu & Wahba 4 term additive model
```

```r
# Rouler un modèle avec intercepte aléatoire
gamm_intercept <- gam(y ~ s(x0) + s(fac, bs = "re"), 
 data = gam_data2, method = "REML")
```

```r
# Voir la sortie du modèle
summary(gamm_intercept)$s.table
#             edf   Ref.df          F     p-value
# s(x0)  2.999782 3.736749   3.608906 0.008655953
# s(fac) 2.977044 3.000000 129.697933 0.000000000
```

---
# GAMM avec un intercepte aléatoire

```r
plot(gamm_intercept, select = 2)
```

---
# GAMM avec un intercepte aléatoire

Nous allons premièrement tracer l'effet combiné de `x0` (sans les niveaux de l'effet aléatoire) et ensuite une courbe pour les 4 niveaux de `fac` :

```r
par(mfrow = c(1,2), cex = 1.1)

# Visualiser les effets combinés de x0 (sans effets aléatoires)
plot_smooth(gamm_intercept, view = "x0", rm.ranef = TRUE,
            main = "intercept + s(x1)")

# Visualiser chaque niveau de l'effet aléatoire
plot_smooth(gamm_intercept, view = "x0", rm.ranef = FALSE,
            cond = list(fac="1"),
            main = "... + s(fac)", col = 'orange', ylim = c(0,25))
plot_smooth(gamm_intercept, view = "x0", rm.ranef = FALSE,
            cond = list(fac = "2"),
            add = TRUE, col = 'red')
plot_smooth(gamm_intercept, view="x0", rm.ranef = FALSE,
            cond = list(fac = "3"),
            add = TRUE, col = 'purple')
plot_smooth(gamm_intercept, view="x0", rm.ranef = FALSE,
            cond = list(fac = "4"),
            add = TRUE, col = 'turquoise')
```

---
# GAMM avec un intercepte aléatoire

---
# GAMM avec une pente aléatoire

Ensuite, spécifions un modèle avec une pente aléatoire:

```r
gamm_slope <- gam(y ~ s(x0) + s(x0, fac, bs = "re"), 
 data = gam_data2, method = "REML")
```

```r
summary(gamm_slope)$s.table
#                edf   Ref.df        F     p-value
# s(x0)     2.942915 3.667881  3.73027 0.008300806
# s(x0,fac) 2.966452 3.000000 88.09058 0.000000000
```

---
# GAMM avec une pente aléatoire

```r
par(mfrow = c(1,2), cex = 1.1)

# Visualiser les effets combinés de x0 (sans effets aléatoires)
plot_smooth(gamm_slope, view = "x0", rm.ranef = TRUE,
            main = "intercept + s(x1)")

# Visualiser chaque niveau de l'effet aléatoire
plot_smooth(gamm_slope, view = "x0", rm.ranef = FALSE,
            cond = list(fac="1"),
            main = "... + s(fac, x0)", col = 'orange', ylim = c(0,25))
plot_smooth(gamm_slope, view = "x0", rm.ranef = FALSE,
            cond = list(fac = "2"),
            add = TRUE, col = 'red')
plot_smooth(gamm_slope, view="x0", rm.ranef = FALSE,
            cond = list(fac = "3"),
            add = TRUE, col = 'purple')
plot_smooth(gamm_slope, view="x0", rm.ranef = FALSE,
            cond = list(fac = "4"),
            add = TRUE, col = 'turquoise')
```

---
# GAMM avec une pente aléatoire

---
# GAMM avec un intercepte et une pente aléatoire

On peut aussi inclure un intercept et une pente aléatoire.

```r
gamm_int_slope <- gam(y ~ s(x0) + s(fac, bs = "re") + s(fac, x0, bs = "re"),
 data = gam_data2, method = "REML")

summary(gamm_int_slope)$s.table
#                edf   Ref.df          F     p-value
# s(x0)     3.006867 3.744708   3.593945 0.008802081
# s(fac)    2.940108 3.000000 381.418397 0.000000000
# s(fac,x0) 1.388805 3.000000 109.266270 0.113010130
```

---
# GAMM avec un intercepte et une pente aléatoire

```r
par(mfrow = c(1,2), cex = 1.1)

# Visualiser les effets combinés de x0 (sans effets aléatoires)
plot_smooth(gamm_int_slope, view = "x0", rm.ranef = TRUE,
            main = "intercept + s(x1)")

# Visualiser chaque niveau de l'effet aléatoire
plot_smooth(gamm_int_slope, view = "x0", rm.ranef = FALSE,
            cond = list(fac="1"),
            main = "... + s(fac) + s(fac, x0)", col = 'orange', ylim = c(0,25))
plot_smooth(gamm_int_slope, view = "x0", rm.ranef = FALSE,
            cond = list(fac = "2"),
            add = TRUE, col = 'red')
plot_smooth(gamm_int_slope, view="x0", rm.ranef = FALSE,
            cond = list(fac = "3"),
            add = TRUE, col = 'purple')
plot_smooth(gamm_int_slope, view="x0", rm.ranef = FALSE,
            cond = list(fac = "4"),
            add = TRUE, col = 'turquoise')
```

---
# GAMM avec un intercepte et une pente aléatoire

---
# GAMM avec un intercepte et une pente aléatoire

Notez que la pente aléatoire est statique dans ce cas :

```r
plot(gamm_int_slope, select = 3)
```

???

select = 3 parce que la pente aléatoire est sur la troisième ligne de notre tableau sommaire.

---
# GAMM avec une surface lisse aléatoire

Finalement, spécifions un modèle avec une surface lisse aléatoire.

```r
gamm_smooth <- gam(y ~ s(x0) + s(x0, 
 fac, 
 bs = "fs",
 m = 1),
 data = gam_data2, 
 method = "REML")
```

```r
summary(gamm_smooth)$s.table
#                edf   Ref.df         F    p-value
# s(x0)     2.973666  3.69337  3.415788 0.01190229
# s(x0,fac) 3.984247 35.00000 11.185386 0.00000000
```

]

```r
plot(gamm_smooth, select=1)
```

]

???
select = 1 parce que la surface lisse aléatoire est sur la première ligne de notre tableau sommaire.

---
# GAMM avec une surface lisse aléatoire

```r
par(mfrow = c(1,2), cex = 1.1)

# Visualiser les effets combinés de x0 (sans effets aléatoires)
plot_smooth(gamm_smooth, view = "x0", rm.ranef = TRUE,
            main = "intercept + s(x1)")

# Visualiser chaque niveau de l'effet aléatoire
plot_smooth(gamm_smooth, view = "x0", 
            rm.ranef = FALSE,
            cond = list(fac="1"),
            main = "... + s(fac) + s(fac, x0)", 
            col = 'orange', ylim = c(0,25))
plot_smooth(gamm_smooth, view = "x0",
            rm.ranef = FALSE,
            cond = list(fac = "2"),
            add = TRUE, col = 'red')
plot_smooth(gamm_smooth, view="x0", 
            rm.ranef = FALSE,
            cond = list(fac = "3"),
            add = TRUE, col = 'purple')
plot_smooth(gamm_smooth, view="x0", 
            rm.ranef = FALSE,
            cond = list(fac = "4"),
            add = TRUE, col = 'turquoise')
```

---
# GAMM avec une surface lisse aléatoire

???

Ici, si la pente aléatoire varie selon `x0`, nous aurions des courbes variables pour chaque niveau.

---
# Comparaison de GAMM

Tous les GAMMs de cette section peuvent être comparé avec `AIC()` pour trouver le modèle le mieux ajusté:

```r
AIC(gamm_intercept, gamm_slope, gamm_int_slope, gamm_smooth)
#                       df      AIC
# gamm_intercept  8.623101 2200.672
# gamm_slope      8.450850 2269.562
# gamm_int_slope 10.840655 2201.538
# gamm_smooth    10.558836 2202.343
```

Le meilleur modèle de ces trois modèles serait donc le GAMM avec un intercept aléatoire.

---
# Ressources

Cet atelier présente une brève introduction aux concepts de base et aux librairies populaires pour vous aider à estimer, évaluer et visualiser les GAMs dans R, mais il y a beaucoup plus à explorer!
 
* Le livre [Generalized Additive Models: An Introduction with R](https://www.routledge.com/Generalized-Additive-Models-An-Introduction-with-R-Second-Edition/Wood/p/book/9781498728331) par Simon Wood (auteur de la librairie `mgcv`).
* Le site web de Simon Wood, [maths.ed.ac.uk/~swood34/](https://www.maths.ed.ac.uk/~swood34/).
* Le blogue de Gavin Simpson, [From the bottom of the heap](https://fromthebottomoftheheap.net/).
* La librairie [`gratia`](https://cran.r-project.org/web/packages/gratia/index.html) de Gavin Simpson for GAM visualisation in `ggplot2`.
* Le cours [Generalized Additive Models: An Introduction with R](https://noamross.github.io/gams-in-r-course/) de Noam Ross.
* [Overview GAMM analysis of time series data](https://jacolienvanrij.com/Tutorials/GAMM.html) tutoriel de Jacolien van Rij.
* [Hierarchical generalized additive models in ecology: an introduction with mgcv](https://peerj.com/articles/6876/) de Pedersen et al. (2019).

Enfin, les pages d'aide, disponibles via `?gam` dans R sont une excellente ressource.

---
class: inverse, center, bottom

# Merci pour votre participation à cet atelier!

<hr>

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

---
class: inverse, center, middle

# Exemple supplémentaire avec d'autres distributions

---
# GAM avec d'autres distributions

Un bref aperçu de l'utilisation des GAMs lorsque la variable réponse ne suit pas une distribution normale ou que les données sont des abondances ou proportions (par exemple, distribution Gamma, binomiale, Poisson, binomiale négative).

Nous allons utiliser un exemple de données où une répartition binomiale sera nécessaire; la variable réponse représente le nombre de succès (l'événement a eu lieu) en fonction des défaillances au cours d'une expérience.

```r
gam_data3 <- read.csv("data/other_dist.csv")
str(gam_data3)
# 'data.frame':	514 obs. of 4 variables:
# $ prop : num 1 1 1 1 0 1 1 1 1 1 ...
# $ total: int 4 20 20 18 18 18 20 20 20 20 ...
# $ x1 : int 550 650 750 850 950 650 750 850 950 550 ...
# $ fac : chr "f1" "f1" "f1" "f1" ...
```

---
# GAM avec d'autres distributions

```r
plot(range(gam_data3$x1), c(0,1), type = "n",
     main = "Probabilités de succès dans le temps",
     ylab = "Probabilité", xlab = "x1 (temps)")
abline(h = 0.5)

avg <- aggregate(prop ~ x1, data=gam_data3, mean)
lines(avg$x1, avg$prop, col = "orange", lwd = 2)
```

---
# GAM avec d'autres distributions

Nous allons tester si cette tendance est linéaire ou non avec un GAM logistique (nous utilisons une famille de distributions binomiales parce que nous avons des données de proportion).

```r
prop_model <- gam(prop ~ s(x1), data = gam_data3, weights = total,
 family = "binomial", method = "REML")
prop_summary <- summary(prop_model)
```

```r
prop_summary$p.table
#             Estimate Std. Error  z value Pr(>|z|)
# (Intercept) 1.174024 0.02709868 43.32402        0
```

```r
prop_summary$s.table
#            edf   Ref.df   Chi.sq p-value
# s(x1) 4.750198 5.791279 799.8463       0
```

---
# GAM avec d'autres distributions

```
#             Estimate Std. Error  z value Pr(>|z|)
# (Intercept) 1.174024 0.02709868 43.32402        0
```

**Rappel** le modèle utilise le nombre de succès vs échecs pour calculer le *logit*, qui est la logarithme du rapport entre les succès et échecs :

.small[
- Si succès = échecs, le rapport = 1 et le logit est de 0 (log(1) = 0).
- Si succès > échecs, le rapport > 1 et le logit a une valeur positive (log(2) = 0.69).
- Si succès < échecs, le rapport < 1 et le logit a une valeur négative (log(.5) = -0.69).
]

> Ici, l'estimé est positif ce qui signifie, qu'en moyenne, il y a plus de succès que d'échecs.

---
# GAM avec d'autres distributions

```
#            edf   Ref.df   Chi.sq p-value
# s(x1) 4.750198 5.791279 799.8463       0
```

Cela représente comment le ratio de succès vs échecs change sur l'échelle de `$x1$`.

> Puisque les **EDF** > 1, la proportion des succès augmente plus rapidement avec `$x1$`

```r
plot(prop_model)
```

---
# Visualiser la tendance au fil du temps

Il y a différente façon de représenter cette relation graphiquement :

- **Contribution/effet partiel** correspond aux effets isolés d'une interaction ou prédiction particulière. Si vous visualisez votre modèle GAM avec `plot()`, vous obtenez les effets partiels.
- **effets additionnés** correspond aux mesures réponse prédites pour une valeur ou niveau donné de prédicteurs. Si vous visualisez votre GAM avec `itsadug::plot_smooth()`, vous obtenez les effets additionnés.

---
# Visualiser la tendance au fil du temps

Que nous disent ces graphes sur les succès et échecs?

La valeur logit augmente, donc les succès augmentent et les échecs diminuent.]

Quantités égales de succès et d'échecs jusqu'à `$x1 = 400$`.
]

---
# Visualiser la tendance au fil du temps

Enfin, pour nous aider à interpréter les résultats, nous pouvons re-transformer l'effet sur une échelle de proportions avec la fonction `itsadug::plot_smooth()` :

```r
plot_smooth(prop_model, view = "x1", main = "",
            transform = plogis, ylim = c(0,1), print.summary = F)
abline(h = 0.5, v = diff$start, col = 'red', lty = 2)
```

Comme précédemment, la proportion de succès augmente au-dessus de 0.5 à `$x1 = 400$`.