11-GAM.qmd

# Modèles généralisés additifs  {#sec-chap11}

Dans les précédents chapitres, nous avons eu l'occasion d'explorer toute une panoplie de modèles : régressions linéaires, modèles généralisés, modèles généralisés à effets mixtes et modèles multiniveaux. Dans ce chapitre, nous abordons une nouvelle extension dans le monde des régressions : les modèles généralisés additifs (*Generalized additive model* en anglais — GAM). Cette extension a pour but de permettre de modéliser des relations non linéaires entre les variables indépendantes et la variable dépendante.

::: bloc_package
::: bloc_package-header
::: bloc_package-icon
:::
**Liste des *packages* utilisés dans ce chapitre**
:::
::: bloc_package-body

* Pour créer des graphiques :
  - `ggplot2` le seul, l'unique!
  - `ggpubr` pour combiner des graphiques et réaliser des diagrammes.
  - `metR` pour placer des étiquettes sur des isolignes.
  
* Pour jouer avec des *splines* :
  - `splines2` pour construire les fonctions de base de nombreuses *splines*.
  - `segmented` pour ajuster des modèles avec des coefficients variant par segment.
  
* Pour ajuster des modèles GAM :
  - `mgcv`, le *package* de référence pour ajuster des GAM dans R!
  - `gamlss`, un second *package* très flexible pour ajuster des GAM.
  - `gamlss.add`, une extension de `gamlss` ajoutant des distributions supplémentaires.

* Pour analyser des modèles GAM :
  - `itsadug` pour notamment extraire certains résultats d'un GAM.
  - `mixedup` pour notamment extraire les effets aléatoires d'un GAM.
  - `DHARMa` pour le diagnostic des résidus simulés.
:::
:::


## Introduction {#sec-111}

Puisque les modèles GAM sont une extension des modèles GLM, ils peuvent s'appliquer à des modèles pour des variables indépendantes qualitatives, de comptage ou continues. Nous l'appliquons ici, à titre d'illustration, à une variable indépendante continue. Pour rappel, la formule décrivant un modèle linéaire généralisé (GLM) utilisant une distribution normale et une fonction de lien identitaire est la suivante : 

$$
\begin{aligned}
&Y \sim Normal(\mu,\sigma)\\
&g(\mu) = \beta_0 + \beta X\\
&g(x) = x
\end{aligned}
$$ {#eq-glm1B}

Les coefficients $\beta$ permettent de quantifier l'effet des variables indépendantes (*X*) sur la moyenne (l'espérance) ($\mu$) de la variable dépendante (*Y*). Un coefficient $\beta_k$ négatif indique que, si la variable $X_k$ augmente, alors la variable *Y* tend à diminuer et inversement, si le coefficient est positif. L'inconvénient de cette formulation est que le modèle est capable de capter uniquement des relations linéaires entre ces variables. Or, il existe de nombreuses situations dans lesquelles une variable indépendante a un lien non linéaire avec une variable dépendante; voici quelques exemples : 

* Si nous mesurons le niveau de bruit émis par une source sonore (variable dépendante) à plusieurs endroits et que nous tentons de prédire l'intensité sonore en fonction de la distance à la source (variable indépendante), nous pouvons nous attendre à observer une relation non linéaire entre les deux. En effet, le son étant une énergie se dispersant selon une sphère dans l'espace, son intensité est inversement proportionnelle au carré de la distance avec la source sonore.

* La concentration de la pollution atmosphérique en ville suit généralement des patrons temporels et spatiaux influencés directement par la météorologie et les activités humaines. Autrement dit, il serait absurde d'introduire l'espace de façon linéaire (avec un gradient nord-sud ou est-ouest), ou le moment de la journée de façon linéaire (comme si la pollution augmentait du matin au soir ou inversement). En guise d'exemple, la @fig-gam1, tirée de @2020_3, illustre bien ces variations temporelles pour deux polluants (le dioxyde d'azote et l'ozone).


![Patron journalier du dioxyde d'azote et de l'ozone à Paris](images/Chap11/no2_03_patterns.png){#fig-gam1 width="75%" fig-align="center"}

### Non-linéarité fonctionnelle {#sec-1111}

Il existe de nombreuses façons d'introduire des relations non linéaires dans un modèle. La première et la plus simple à mettre en œuvre est de transformer la variable indépendante à l'aide d'une fonction inverse, exponentielle, logarithmique ou autre.

Prenons un premier exemple avec une variable *Y* que nous tentons de prédire avec une variable *X*, présenté à la @fig-gam2. Si nous ajustons une droite de régression à ces données (en bleu), nous constatons que l'augmentation de *X* est associée à une augmentation de *Y*. Cependant, la droite de régression est très éloignée des données et ne capte qu'une petite partie de la relation. Une lecture attentive permet de constater que l'effet de *X* sur *Y* augmente de plus en plus rapidement à mesure que *X* augmente. Cette forme est caractéristique d'une relation exponentielle. Nous pouvons donc transformer la variable *X* avec la fonction exponentielle afin d'obtenir un meilleur ajustement (en rouge).

```{r}
#| label: fig-gam2
#| fig-cap: Relation non linéaire exponentielle
#| echo: false
#| fig-align: center
#| message: false
#| warning: false
#| out-width: "50%"
library(ggplot2)
library(ggpubr)
library(boot)

tofr <- function(float){
  return(gsub("." , ",", as.character(float), fixed = TRUE,  useBytes = TRUE))
}

df <- data.frame(
  x1 = rnorm(1000,0,1),
  x2 = rnorm(1000,10,5),
  x3 = rnorm(1000,0,3),
  x4 = abs(rnorm(1000,0,5))
)

df$y1 <- exp(df$x1) + rnorm(1000,0,2)
df$y2 <- -log(df$x2 + rnorm(1000,0,1))
df$y3 <- inv.logit(df$x3) + rnorm(1000,0,0.3)
df$y4 <- sqrt(df$x4) + rnorm(1000,0,0.2)


ggplot(df)+
  geom_point(aes(x = x1, y = y1), size = 0.5) + 
  labs(x = "x", y = "y")  + 
  geom_smooth(aes(x = x1, y = y1), method = "lm",
              se = FALSE, color = "blue", formula = y ~ x)+ 
  geom_smooth(aes(x = x1, y = y1), method = "lm",
              se = FALSE, color = "red", formula = y ~ exp(x))

```


La @fig-gam3 illustre trois autres situations avec les fonctions logarithmique, logistique inverse et racine carrée. Cette approche peut donner des résultats intéressants si vous disposez d'une bonne justification théorique sur la forme attendue de la relation entre *X* et *Y*. 

```{r}
#| label: fig-gam3
#| fig-cap: Autres relations non linéaires
#| echo: false
#| fig-align: center
#| message: false
#| warning: false
#| out-width: "75%"


P1 <- ggplot(df)+
  geom_point(aes(x = x2, y = y2), size = 0.5) + 
  labs(x = "x", y = "y", subtitle = "fonction logarithmique")  + 
  geom_smooth(aes(x = x2, y = y2), method = "lm",
              se = FALSE, color = "blue", formula = y ~ x)+ 
  geom_smooth(aes(x = x2, y = y2), method = "lm",
              se = FALSE, color = "red", formula = y ~ log(x)) + 
  ylim(-3,2.5) + xlim(0,20)

P2 <- ggplot(df)+
  geom_point(aes(x = x3, y = y3), size = 0.5) + 
  labs(x = "x", y = "y", subtitle = "fonction logistique inverse")  + 
  geom_smooth(aes(x = x3, y = y3), method = "lm",
              se = FALSE, color = "blue", formula = y ~ x)+ 
  geom_smooth(aes(x = x3, y = y3), method = "lm",
              se = FALSE, color = "red", formula = y ~ inv.logit(x))

P3 <- ggplot(df)+
  geom_point(aes(x = x4, y = y4), size = 0.5) + 
  labs(x = "x", y = "y", subtitle = "fonction racine carrée")  + 
  geom_smooth(aes(x = x4, y = y4), method = "lm",
              se = FALSE, color = "blue", formula = y ~ x)+ 
  geom_smooth(aes(x = x4, y = y4), method = "lm",
              se = FALSE, color = "red", formula = y ~ sqrt(x))

ggarrange(P1, P2, P3, ncol = 2, nrow = 2)

```

Il existe également des cas de figure dans lesquels aucune fonction ne donne de résultats pertinents, comme illustré à la @fig-gam4. Nous constatons facilement qu'aucune des fonctions proposées n'est capable de bien capter la relation entre les deux variables. Puisque cette relation est complexe, il convient alors d'utiliser une autre stratégie pour la modéliser. 

```{r}
#| label: fig-gam4
#| echo: false
#| fig-align: center
#| fig-cap: Relation non linéaire plus complexe
#| message: false
#| warning: false
#| out-width: "75%"

library(mgcv)
dataset <- gamSim(eg = 1, n = 400, dist = "normal", scale = 1, verbose = FALSE)

ggplot(dataset) + 
  geom_point(aes(y = y, x = x2), size = 1) + 
  geom_smooth(aes(x = x2, y = y, color = "lineaire"), method = "lm",
              se = FALSE, formula = y ~ x, size = 1.3) + 
  geom_smooth(aes(x = x2, y = y, color = "logarithme"), method = "lm",
              se = FALSE, formula = y ~ log(x), size = 1.3) +
  geom_smooth(aes(x = x2, y = y, color = "racine carrée"), method = "lm",
              se = FALSE, formula = y ~ sqrt(x), size = 1.3) + 
  scale_color_manual(values = c("lineaire" = "#0077b6",
                                "logarithme" = "#e63946",
                                "racine carrée" = "#2a9d8f"))+
  labs(colour = 'forme fonctionelle')

```


### Non-linéarité avec des polynomiales {#sec-1112}

Nous avons vu, dans le chapitre sur la régression simple ([section @sec-07511]), qu'il est possible d'utiliser des polynomiales pour ajuster des relations non linéaires. Pour rappel, il s'agit simplement d'ajouter à un modèle la variable *X* à différents exposants ($X+X^2+\dots+X^k$). Chaque exposant supplémentaire (chaque ordre supplémentaire) permet au modèle d'ajuster une relation plus complexe. Rien de tel qu'un graphique pour illustrer le tout (@fig-gam5).

```{r}
#| label: fig-gam5
#| echo: false
#| fig-align: center
#| fig-cap: Visualisation de plusieurs polynomiales
#| message: false
#| warning: false
#| out-width: "70%"

p1 <- ggplot(dataset) + 
  geom_point(aes(y = y, x = x2), size = 1) + 
  geom_smooth(aes(x = x2, y = y), color = "blue" , method = "lm",
              se = FALSE, formula = y ~ x + I(x**2), size = 1.3) + 
  labs(subtitle = "polynomiale de degré 2", x = "x", y = "y")
  
p2 <- ggplot(dataset) + 
  geom_point(aes(y = y, x = x2), size = 1) + 
  geom_smooth(aes(x = x2, y = y), color = "blue" , method = "lm",
              se = FALSE, formula = y ~ x + I(x**2) + I(x**3), size = 1.3) + 
  labs(subtitle = "polynomiale de degré 3", x = "x", y = "y")

p3 <- ggplot(dataset) + 
  geom_point(aes(y = y, x = x2), size = 1) + 
  geom_smooth(aes(x = x2, y = y), color = "blue" , method = "lm",
              se = FALSE, formula = y ~ x + I(x**2) + I(x**3) + I(x**4), size = 1.3) + 
  labs(subtitle = "polynomiale de degré 4", x = "x", y = "y")

p4 <- ggplot(dataset) + 
  geom_point(aes(y = y, x = x2), size = 1) + 
  geom_smooth(aes(x = x2, y = y), color = "blue" , method = "lm",
              se = FALSE, formula = y ~ x + I(x**2) + I(x**3) + I(x**4)+ I(x**5), size = 1.3) + 
  labs(subtitle = "polynomiale de degré 5", x = "x", y = "y")

ggarrange(p1, p2, p3, p4, ncol = 2, nrow = 2)
```

L'enjeu est de sélectionner le bon nombre de degrés de la polynomiale pour le modèle. Chaque degré supplémentaire constitue une nouvelle variable dans le modèle, et donc un paramètre supplémentaire. Un trop faible nombre de degrés produit des courbes trop simplistes, alors qu'un nombre trop élevé conduit à un surajustement (*overfitting* en anglais) du modèle. La @fig-gam6 illustre ces deux situations.

```{r}
#| label: fig-gam6
#| echo: false
#| fig-align: center
#| fig-cap: Sur et sous-ajustement d'une polynomiale
#| message: false
#| warning: false
#| out-width: "70%"

ggplot(dataset) + 
  geom_point(aes(y = y, x = x2), size = 1) + 
  geom_smooth(aes(x = x2, y = y), color = "blue" , method = "lm",
              se = FALSE, formula = y ~ x + I(x**2), size = 1.3) +
  geom_smooth(aes(x = x2, y = y), color = "red" , method = "lm",
              se = FALSE, formula = y ~ x + poly(x, degree = 12), size = 1.3) + 
  labs(subtitle = "polynomiales de degrés 2 (bleu) et 12 (rouge)", x = "x", y = "y")
```

Un des problèmes inhérents à l'approche des polynomiales est la difficulté d'interprétation. En effet, les coefficients ne sont pas directement interprétables et seule une figure représentant les prédictions du modèle permet d'avoir une idée de l'effet de la variable *X* sur la variable *Y*.

### Non-linéarité par segments {#sec-1113}

Un compromis intéressant offrant une interprétation simple et une relation potentiellement complexe consiste à découper la variable *X* en segments, puis d'ajuster un coefficient pour chacun de ces segments. Nous obtenons ainsi une ligne brisée et des coefficients faciles à interpréter (@fig-gam7). Nous ne présentons pas d'exemple d'application dans R, mais sachez que le *package* `segmented` permet d'ajuster ce type de modèle. 

```{r}
#| label: fig-gam7
#| echo: false
#| fig-align: center
#| fig-cap: Régression par segment
#| message: false
#| warning: false
#| out-width: "70%"

library(segmented)

model <- lm(y ~ x2, data = dataset)
o <- segmented(model, seg.Z = ~x2, psi = list(x2 = c(0.25,0.5)),
  control = seg.control(display = FALSE)
)
dataset$fit <- broken.line(o)$fit

ggplot(dataset) + 
  geom_point(aes(y = y, x = x2), size = 1) + 
  geom_line(aes(x = x2, y = fit), color = "blue", linewidth = 1)
```

L'enjeu est alors de déterminer le nombre de points et la localisation de points de rupture. L'inconvénient majeur de cette approche est qu'en réalité, peu de phénomènes sont marqués par des ruptures très nettes.

À la @fig-gam7, nous avons divisé la variable *X* en trois segments ($k_1$, $k_2$ et $k_3$), définis respectivement avec les intervalles suivants : [0,00-0,22], [0,22-0,41] et [0,41-1,00]. Concrètement, cela revient à diviser la variable *X* en trois nouvelles variables $X_{k1}$, $X_{k2}$, et $X_{k3}$. La valeur de $X_{ik}$ est égale à $x_i$ si $x_i$ se trouve dans l'intervalle propre à *k*, et à 0 autrement. Ici, nous obtenons trois coefficients : 

* le premier est positif, une augmentation de *X* sur le premier segment est associée à une augmentation de *Y*;

* le second est négatif, une augmentation de *X* sur le second segment est associée à une diminution de *Y*;

* le troisième est aussi négatif, une augmentation de *X* sur le troisième segment est associée à une diminution de *Y*, mais moins forte que pour le second segment.


### Non-linéarité avec des *splines* {#sec-1114}

La dernière approche, et certainement la plus flexible, est d'utiliser ce que l'on appelle une *spline* pour capter des relations non linéaires. Une *spline* est une fonction créant des variables supplémentaires à partir d'une variable *X* et d'une fonction de base. Ces variables supplémentaires, appelées bases (*basis* en anglais), sont ajoutées au modèle; la sommation de leurs valeurs multipliées par leurs coefficients permet de capter les relations non linéaires entre une variable dépendante et une variable indépendante. Le nombre de bases et leur localisation (plus souvent appelé nœuds) permettent de contrôler la complexité de la fonction non linéaire.

Prenons un premier exemple simple avec une fonction de base triangulaire (*tent basis* en anglais). Nous créons ici une *spline* avec sept nœuds répartis équitablement sur l'intervalle de valeurs de la variable *X*. Les sept bases qui en résultent sont présentées à la @fig-gam8. Dans cette figure, chaque sommet d'un triangle correspond à un nœud et chaque triangle correspond à une base.


```{r}
#| label: fig-gam8
#| echo: false
#| fig-align: center
#| fig-cap: Bases de la spline triangulaire
#| message: false
#| warning: false
#| out-width: "70%"
library(splines)

basis <- bs(dataset$x2, df = 7, degre = 1)
df <- data.frame(basis * 3)
df$X <- dataset$x2
df <- reshape2::melt(df, id.vars = "X")

ggplot() + 
  geom_point(aes(y = y, x = x2), size = 1, data = dataset) + 
  geom_line(aes(x = X, y = value, color = variable), linewidth = 1, data = df) + theme(legend.position = "none") 
```


En ajoutant ces bases dans notre modèle de régression, nous pouvons ajuster un coefficient pour chacune et le représenter en multipliant ces bases par les coefficients obtenus avec une simple régression linéaire (@fig-gam9).

```{r}
#| label: fig-gam9
#| echo: false
#| fig-align: center
#| fig-cap: Spline triangulaire multipliée par ces coefficients
#| message: false
#| warning: false
#| out-width: "70%"

basis <- bs(dataset$x2, df = 7, degre = 1)
model <- lm(y ~ basis, data = dataset)
coefs <- coefficients(model)[2:length(coefficients(model))]

prodbase <- sweep(basis, MARGIN = 2, coefs, `*`)

df <- data.frame(prodbase)
df$X <- dataset$x2
df <- reshape2::melt(df, id.vars = "X")

ggplot() + 
  geom_point(aes(y = y, x = x2), size = 1, data = dataset) + 
  geom_line(aes(x = X, y = value, color = variable), linewidth = 1, data = df) + theme(legend.position = "none") 
```

Nous remarquons ainsi que les bases correspondant à des valeurs plus fortes de *Y* ont reçu des coefficients plus élevés. Pour reconstituer la fonction non linéaire, il suffit d'additionner ces bases multipliées par leurs coefficients, soit la ligne bleue à la @fig-gam10.

```{r}
#| label: fig-gam10
#| echo: false
#| fig-align: center
#| fig-cap: Spline triangulaire
#| message: false
#| warning: false
#| out-width: "70%"


basis <- bs(dataset$x2, df = 7, degre = 1)
model <- lm(y ~ basis, data = dataset)
coefs <- coefficients(model)[2:length(coefficients(model))]

prodbase <- sweep(basis, MARGIN = 2, coefs, `*`)
total <- rowSums(prodbase) + coefficients(model)[1]

df <- data.frame(prodbase)
df$X <- dataset$x2
df <- reshape2::melt(df, id.vars = "X")
df$total <- total

ggplot() + 
  geom_point(aes(y = y, x = x2), size = 1, data = dataset) + 
  geom_line(aes(x = X, y = value, group = variable), color = "grey", linetype = "dashed", linewidth = 1, data = df) +
  geom_line(aes(x = X, y = total), color = "blue", linewidth = 1, data = df)+
  theme(legend.position = "none") 
```

La fonction de base triangulaire est intéressante pour présenter la logique qui sous-tend les *splines*, mais elle est rarement utilisée en pratique. On lui préfère généralement d'autres formes donnant des résultats plus lisses comme les *B-spline* quadratiques, *B-spline* cubiques, *M-spline*, *Duchon spline*, etc.

```{r}
#| label: fig-gam11
#| echo: false
#| fig-align: center
#| fig-cap: Comparaison de différentes bases
#| message: false
#| warning: false
#| out-width: "75%"

library(mgcv)
library(splines2)

# Squared spline
basis <- bSpline(dataset$x2, df = 7, degre = 2)
model <- lm(y ~ basis, data = dataset)
coefs <- coefficients(model)[2:length(coefficients(model))]

prodbase <- sweep(basis, MARGIN = 2, coefs, `*`)
total <- rowSums(prodbase) + coefficients(model)[1]

df <- data.frame(prodbase)
df$X <- dataset$x2
df <- reshape2::melt(df, id.vars = "X")
df$total <- total

P1 <- ggplot() + 
  geom_point(aes(y = y, x = x2), size = 1, data = dataset) + 
  geom_line(aes(x = X, y = value, group = variable), color = "grey", linetype = "dashed", linewidth = 1, data = df) +
  geom_line(aes(x = X, y = total), color = "blue", linewidth = 1, data = df)+
  theme(legend.position = "none") + 
  labs(subtitle = "B-spline quadratique")+ 
  ylim(0,20)


# Spline cubique
basis <- bSpline(dataset$x2, df = 7, degre = 3)
model <- lm(y ~ basis, data = dataset)
coefs <- coefficients(model)[2:length(coefficients(model))]

prodbase <- sweep(basis, MARGIN = 2, coefs, `*`)
total <- rowSums(prodbase) + coefficients(model)[1]

df <- data.frame(prodbase)
df$X <- dataset$x2
df <- reshape2::melt(df, id.vars = "X")
df$total <- total

P2 <- ggplot() + 
  geom_point(aes(y = y, x = x2), linewidth = 1, data = dataset) + 
  geom_line(aes(x = X, y = value, group = variable), color = "grey", linetype = "dashed", linewidth = 1, data = df) +
  geom_line(aes(x = X, y = total), color = "blue", linewidth = 1, data = df)+
  theme(legend.position = "none") + 
  labs(subtitle = "B-spline cubique")+ 
  ylim(0,20)

# M-spline
basis <- mSpline(dataset$x2, df = 7, degre = 2)
model <- lm(y ~ basis, data = dataset)
coefs <- coefficients(model)[2:length(coefficients(model))]

prodbase <- sweep(basis, MARGIN = 2, coefs, `*`)
total <- rowSums(prodbase) + coefficients(model)[1]

df <- data.frame(prodbase)
df$X <- dataset$x2
df <- reshape2::melt(df, id.vars = "X")
df$total <- total

P3 <- ggplot() + 
  geom_point(aes(y = y, x = x2), size = 1, data = dataset) + 
  geom_line(aes(x = X, y = value, group = variable), color = "grey", linetype = "dashed", linewidth = 1, data = df) +
  geom_line(aes(x = X, y = total), color = "blue", linewidth = 1, data = df)+
  theme(legend.position = "none") + 
  labs(subtitle = "M-spline")+ 
  ylim(0,20)


# Duchon spline
smoother <- smoothCon(s(x2, bs="ds", k = 7), data = dataset, absorb.cons = T)
basis <- smoother[[1]]$X
model <- lm(y ~ basis, data = dataset)
coefs <- coefficients(model)[2:length(coefficients(model))]

prodbase <- sweep(basis, MARGIN = 2, coefs, `*`)
total <- rowSums(prodbase) + coefficients(model)[1]

df <- data.frame(prodbase)
df$X <- dataset$x2
df <- reshape2::melt(df, id.vars = "X")
df$total <- total

P4 <- ggplot() + 
  geom_point(aes(y = y, x = x2), size = 1, data = dataset) + 
  geom_line(aes(x = X, y = value, group = variable), color = "grey", linetype = "dashed", linewidth = 1, data = df) +
  geom_line(aes(x = X, y = total), color = "blue", linewidth = 1, data = df)+
  theme(legend.position = "none") + 
  labs(subtitle = "Duchon-spline") + 
  ylim(0,20)

ggarrange(P1, P2, P3, P4, ncol = 2, nrow = 2)
```

Les approches que nous venons de décrire sont regroupées sous l'appellation de modèles additifs. Dans les prochaines sous-sections, nous nous concentrons davantage sur les *splines* du fait de leur plus grande flexibilité.

## *Spline* de régression et *spline* de lissage {#sec-112}

Dans les exemples précédents, nous avons vu que la construction d'une *spline* nécessite d'effectuer deux choix importants : le nombre de nœuds et leur localisation. Un trop grand nombre de nœuds conduit à un surajustement du modèle alors qu'un trop faible nombre de nœuds conduit à un sous-ajustement. Lorsque ces choix sont effectués par l'utilisateur et que les bases sont ajoutées manuellement dans le modèle tel que décrit précédemment, nous parlons alors de **_splines_ de régression** (*Regression Spline* en anglais).

Une approche a été proposée pour faciliter le choix du nombre de nœuds, il s'agit de **_splines_ de lissage**  (*smoothing spline* en anglais). L'idée derrière cette approche est d'introduire dans le modèle une pénalisation associée avec le nombre de nœuds (ou degré de liberté) de la *spline*, dans un souci de parcimonie : chaque nœud supplémentaire doit suffisamment contribuer au modèle pour être conservé. Il n'est pas nécessaire ici de rentrer dans le détail mathématique de cette pénalisation qui est un peu complexe. Retenez simplement qu'elle dépend d'un paramètre appelé $\lambda$ :

* plus $\lambda$ tend vers 0, plus la pénalisation est faible et plus la *spline* de lissage devient une simple *spline* de régression; 

* à l'inverse, plus elle est forte, plus la pénalité est importante, au point que la *spline* peut se résumer à une simple ligne droite. 

Cela est illustré à la @fig-gam12 comprenant trois *splines* avec 20 nœuds et des valeurs $\lambda$ différentes contrôlant la force de la pénalité. 

Bien évidemment, nous constatons qu'avec la *spline* de régression (non pénalisée), 20 nœuds conduisent à un fort surajustement du modèle. En revanche, les *splines* de lissage (pénalisées) permettent de corriger ce problème de surajustement. Toutefois, une valeur trop importante de $\lambda$ conduit à un sous-ajustement du modèle (ici $\lambda = 3$ et $\lambda = 100$, lignes verte et bleue).


```{r}
#| label: fig-gam12
#| echo: false
#| fig-align: center
#| fig-cap: Pénalisation des splines
#| message: false
#| warning: false
#| out-width: "70%"

# Modèle 1 : régression spline sans pénalite
mod1 <- gam(y~ s(x2, k = 20, fx = T), data = dataset)

# Modèle 2 : régression spline avec pénalite = 3
mod2 <- gam(y~ s(x2, k = 20, sp = 0.5), data = dataset)

# Modèle 3 : régression spline avec pénalite = 10
mod3 <- gam(y~ s(x2, k = 20, sp = 3), data = dataset)

# Modèle 4 : régression spline avec pénalite = 100
mod4 <- gam(y~ s(x2, k = 20, sp = 100), data = dataset)

df <- data.frame(
  y = dataset$y,
  x = dataset$x2,
  pred1 = predict(mod1),
  pred2 = predict(mod2),
  pred3 = predict(mod3),
  pred4 = predict(mod4)
)

ggplot(df) + 
  geom_point(aes(x = x, y = y), size = 1) +
  geom_line(aes(x = x, y = pred1, color = "spline de régression"), linewidth = 1.1) +
  geom_line(aes(x = x, y = pred2, color = "lambda = 0.5"), linewidth = 1.1) +
  geom_line(aes(x = x, y = pred3, color = "lambda = 3"), linewidth = 1.1)+
  geom_line(aes(x = x, y = pred4, color = "lambda = 100"), linewidth = 1.1) + 
  labs(colour = 'pénalisation')

```


Avec les *splines* de lissage, l'enjeu est de sélectionner une valeur optimale de $\lambda$. Le plus souvent, les *packages* R **estiment eux-mêmes** ce paramètre à partir des données utilisées dans le modèle. Toutefois, gardez en mémoire que vous pouvez modifier ce paramètre. Mentionnons également que les *splines* de lissage peuvent être reparamétrées dans un modèle pour être intégrées comme des effets aléatoires. Dans ce cas-ci, $\lambda$ est remplacé par un simple paramètre de variance directement estimé dans le modèle [@wood2004stable].


## Interprétation d'une *spline* {#sec-113}

L'interprétation d'une *spline* se fait à l'aide de graphiques. En effet, puisqu'elle est composée d'un ensemble de coefficients appliqués à des bases, il est difficile d'interpréter directement ces derniers. Nous préférons alors représenter la fonction obtenue à l'aide d'un graphique, illustrant son **effet marginal**. Ce graphique est construit en trois étapes : 

1. Créer un jeu de données fictif dans lequel l'ensemble des variables indépendantes sont fixées à leurs moyennes respectives, sauf la variable pour laquelle nous souhaitons représenter la *spline*. Pour cette dernière, un ensemble de valeurs allant de son minimum à son maximum est utilisé;

2. Utiliser le modèle pour prédire les valeurs attendues de la variable dépendante pour chacune des observations fictives ainsi créées;

3. Afficher les prédictions obtenues dans un graphique.

Notez ici qu'un graphique des effets marginaux se base sur les prédictions du modèle. Si un modèle est mal ajusté, les prédictions ne seront pas fiables et il sera inutile d'interpréter la *spline* obtenue.

Il est aussi possible, dans le cas des *splines* de lissage, d'interpréter les *estimated degrees of freedom* (EDF) qui constituent une approximation du nombre de nœuds de la *spline*. S'ils ne nous renseignent pas sur la forme de la *spline*, ils nous indiquent son niveau de complexité. Une *spline* avec un EDF de 1 est en réalité un simple terme linéaire. Plus l'EDF augmente, plus la *spline* est complexe. 


## Multicolinéarité non linéaire {#sec-114}

Lorsque des *splines* sont ajoutées dans un modèle, il est nécessaire de vérifier si ces dernières ne posent pas un problème de multicolinéarité. Cependant, le VIF ne peut plus être utilisé du fait de la non-linéarité des relations modélisées. Il est alors nécessaire d'utiliser une autre mesure : la concurvité (*concurvity*) permettant de mesurer sur une échelle allant de 0 à 1 à quel point deux *splines* ont en réalité capturé le même effet et se substituent l'une à l'autre. Une valeur de 0 indique une absence totale de concurvité alors qu'une valeur de 1 indique que deux *splines* sont rigoureusement identiques (modèle non identifiable). 

## *Splines* avancées {#sec-115}

Jusqu'ici, nous avons seulement présenté le cas le plus simple pour lequel une *spline* est construite à partir d'une seule variable dépendante continue, mais les **splines** peuvent être utilisées dans de nombreux autres contextes et ont une incroyable flexibilité. Nous détaillons ici trois exemples fréquents : les *splines* cycliques, les *splines* variant par groupe et les *splines* multivariées. Pour une description complète des effets non linéaires possibles avec `mgcv`, n'hésitez pas à consulter sa [documentation](https://stat.ethz.ch/R-manual/R-patched/library/mgcv/html/smooth.terms.html).

```{r}
#| label: tbl-exemplesplines
#| tbl-cap: Exemples de splines avancées
#| echo: false
#| message: false
#| warning: false

df <- data.frame(
        Type = c("spline cyclique", "spline variant par groupe", "spline bivariée", 
                 "spline d'interaction complète", "spline d'interaction partielle"),
        Code = c("`s(x, bs = 'cc')`" , "`s(x, by = x2)`" , "`s(x1,x2)`" , "`te(x1,x2)`",
                 "`s(x1) + s(x2) + ti(x1,x2)`"),
        Description = c("Une spline cyclique doit être utilisée si le 0 de la variable X correspond également à sa valeur maximum. Un bon exemple est le temps dans une journée, car 24 h est équivalent à 0 h", "Une spline variant par groupe permet d'ajuster une spline à une variable X1 différente pour chaque groupe identifié par une variable qualitative X2", "Une spline bivariée est utilisée pour modéliser l'interaction non linéaire de deux variables X1 et X2 s'exprimant dans la même unité (typiquement des coordonnées géographiques cartésiennes)", "Une spline d'interaction permet de modéliser l'interaction non linéaire de deux variables continues pouvant s'exprimer dans des unités différentes, elle combine les effets spécifiques de chacune des deux variables et leur interaction", "Une spline d'interaction partielle permet de distinguer les effets non linéaires individuels de deux variables de leur interaction non linéaire"))

options(knitr.kable.NA = "")
knitr::kable(df, 
           format.args = list(decimal.mark = ',', big.mark = " "),
           col.names = c("Type" , "Code" , "Description"),
           col.to.resize = 3,
           col.width = c("10cm")
           )

```


### *Splines* cycliques {#sec-1151}

Une *spline* cyclique est une extension d'une *spline* classique dont les bases aux extrémités sont spécifiées de telle sorte que la valeur au départ de la *spline* soit la même que celle à la fin de la *spline*. Cela permet à la *spline* de former une boucle, ce qui est particulièrement intéressant pour des variables dont le 0 et la valeur maximale correspondent en réalité à la même valeur. L'exemple le plus parlant est certainement le cas d'une variable représentant la mesure d'un angle en degrés. Les valeurs de 0 et 360 sont identiques et les valeurs 350 et 10 sont toutes les deux à une distance de 10 degrés de 0. Un autre exemple possible serait de considérer l'heure comme une variable continue; dans ce cas, 24 h et 0 h signifient la même chose.

Prenons un exemple concret. Nous souhaitons modéliser la concentration de dioxyde d'azote (NO~2~) à Paris, mesurée par un ensemble de stations fixes. Nous pourrions nous attendre à ce que le NO~2~ suive chaque jour un certain patron. Concrètement, à proximité d'axes routiers majeurs, nous nous attendons à observer des pics suivant les flux pendulaires. À la @fig-gam13, nous retrouvons bien les deux pics attendus correspondant aux heures de pointe du matin et du soir. Aussi, comme indiqué par la ligne rouge, la valeur prédite par la *spline* est la même à 24 h et à 0 h.


```{r}
#| label: fig-gam13
#| echo: false
#| fig-align: center
#| fig-cap: Spline cyclique pour modéliser la concentration de dioxyde d'azote
#| message: false
#| warning: false
#| out-width: "70%"

dataset <- read.csv('data/gam/NO2_airparif.csv', sep=';')
dataset <- dataset[2:nrow(dataset),]
df <- reshape2::melt(dataset, id.vars = c("heure" , "date"))
df$no2 <- ifelse(df$value == "n/d", NA, as.numeric(df$value))

model <- gam(no2 ~ s(heure, bs = "cp"), data = df, family = tw)
preds <- data.frame(
  heure = 1:24
)
modout <- predict(model, newdata = preds, type = "link", se.fit = TRUE)
preds$yhat <- exp(modout$fit)
preds$lower <- exp(modout$fit - 1.96*modout$se.fit)
preds$upper <- exp(modout$fit + 1.96*modout$se.fit)

subpts <- df[sample(1:nrow(df), size = 1000, replace = FALSE),]
subpts$heure2 <- subpts$heure + runif(nrow(subpts), min = -0.5, max = +0.5)

cent_line <- preds$yhat[[1]]

ggplot() + 
  geom_point(mapping = aes(x = heure2, y = no2), size = 0.5, data = subpts, alpha = 0.5) +
  geom_ribbon(mapping = aes(x = heure, ymax = upper, ymin = lower), fill = "grey", alpha = 0.4, data = preds) +
  geom_line(mapping = aes(x = heure, y = yhat), color = "blue", linewidth = 1.1, data = preds) + 
  geom_hline(yintercept = cent_line, color = "red", linetype = "dashed") +
  labs(x = "moment de la journée", y = "concentration de NO2") + 
  scale_x_continuous(breaks = seq(4,20,4), labels = paste(seq(4,20,4),'h', sep='')) + 
  xlim(1,24)

```



### Splines par groupe {#sec-1152}

Tel qu'abordé dans les chapitres précédents, il arrive régulièrement que les observations appartiennent à  différents groupes. Dans ce cas de figure, nous pouvons être amenés à vérifier si la relation décrite par une *spline* est identique pour chacun des groupes d'observations. Il s'agit alors d'ajuster une *spline* différente par groupe. Dans l'exemple précédent, chaque valeur de NO~2~ a été mesurée par une station fixe de mesure spécifique. Compte tenu du fait que l'environnement autour de chaque station est particulier, nous pourrions s'attendre à ce que les valeurs de NO~2~ ne présentent pas exactement les mêmes patrons journaliers pour chaque station.

À la @fig-gam14, il est possible de constater que le NO~2~ suit globalement le même patron temporel pour l'ensemble des stations à l'exception de trois d'entres-elles. Il s'agit en réalité de stations situées dans des secteurs ruraux de la région parisienne, et donc moins impactées par le trafic routier.


```{r}
#| label: fig-gam14
#| echo: false
#| fig-align: center
#| fig-cap: Spline cyclique variant par groupe
#| message: false
#| warning: false
#| out-width: "70%"

model <- gam(no2 ~ s(heure, bs = "cp", by = variable), data = df, family = gaussian)
preds <- expand.grid(
  heure = 1:24,
  variable = unique(df$variable)
)

modout <- predict(model, newdata = preds, type = "link", se.fit = TRUE)
preds$yhat <- modout$fit
preds$lower <- modout$fit - 1.96*modout$se.fit
preds$upper <-modout$fit + 1.96*modout$se.fit

subpts <- df[sample(1:nrow(df), size = 1000, replace = FALSE),]
subpts$heure2 <- subpts$heure + runif(nrow(subpts), min = -0.5, max = +0.5)

cent_line <- preds$yhat[[1]]

ggplot() + 
  geom_point(mapping = aes(x = heure2, y = no2), size = 0.5, data = subpts, alpha = 0.5) +
  geom_ribbon(mapping = aes(x = heure, ymax = upper, ymin = lower, group = variable), fill = "grey", alpha = 0.4, data = preds) +
  geom_line(mapping = aes(x = heure, y = yhat, color = variable), size =0.7, data = preds) + 
  labs(x = "moment de la journée", y = "concentration de NO2", color = "station de mesure") + 
  scale_x_continuous(breaks = seq(4,20,4), labels = paste(seq(4,20,4),'h', sep='')) + 
  xlim(1,24) + 
  theme(legend.position = 'NONE')

```


### *Splines* multivariées et *splines* d'interaction {#sec-1153}

Jusqu'ici, nous n'avons considéré que des *splines* ne s'appliquant qu'à une seule variable indépendante; cependant, il est possible de construire des *splines* multivariées s'ajustant simultanément sur plusieurs variables indépendantes. L'objectif est alors de modéliser les potentielles interactions non linéaires entre les variables indépendantes combinées dans une même *spline*. Prenons un exemple concret, dans la section sur les modèles GLM, nous avons modélisé la couverture des aires de diffusion (AD) à Montréal par des îlots de chaleur. Parmi les variables indépendantes, nous avons notamment utilisé la distance au centre-ville ainsi que la part de la surface végétalisée des AD. Nous pourrions formuler l'hypothèse que ces deux variables influencent conjointement et de façon non linéaire la proportion de la surface d'îlot de chaleur dans chaque AD. Pour représenter une *spline* sur plusieurs dimensions, nous utilisons alors une carte de chaleur dont la couleur représente la valeur de la variable dépendante prédite en fonction des deux variables indépendantes.

Il est important de distinguer la *spline* d'**interaction** et la *spline* **multivariée**. La première est utilisée lorsque les variables indépendantes introduites dans la *spline* ne sont pas exprimées sur la même échelle et n'évoluent pas conjointement. L'exemple donné ci-dessus avec les variables de végétation et de distance au centre-ville est un exemple de *spline* d'interaction, la première variable étant exprimée en pourcentage et l'autre en mètres. De plus, ces deux variables ne sont pas conjointes, mais bien distinctes l'une de l'autre. Un cas typique où une *spline* multivariée serait à privilégier est le cas de l'ajout des coordonnées spatiales dans le modèle. L'emplacement des AD est mesuré par deux variables (coordonnées spatiales *x* et *y*) toutes les deux exprimées en mètres évoluant conjointement, au sens où les coordonnées *x* n'interagissent pas avec les coordonnées *y*, mais forment à elles deux un espace propre. Au-delà de la problématique de l'échelle des données, il est important de retenir que les *splines* d'interaction tendent à être davantage pénalisées que les splines *multivariées*.

La *spline* d'interaction représentée à la @fig-gam15 indique que les AD avec la plus grande proportion de leur surface couverte par des îlots de chaleur sont situées à moins de 25 kilomètres du centre-ville, au-delà de cette distance, cette proportion chute en bas de 0,1, soit 10 % de la surface de l'AD. En revanche, à proximité du centre-ville (moins d'un kilomètre), même les AD disposant d'un fort pourcentage de surface végétalisée sont tout de même marquées par un fort pourcentage de surface couverte par des îlots de chaleur.

Les *splines* bivariées sont fréquemment utilisées pour capturer un potentiel patron spatial dans les données. En effet, si nous disposons des coordonnées spatiales de chaque observation (*x, y*), il est possible d'ajuster une *spline* bivariée sur ces coordonnées, contrôlant ainsi l'effet de l'espace. 

```{r}
#| label: fig-gam15
#| echo: false
#| fig-align: center
#| fig-cap: Spline d'interaction bivariée
#| message: false
#| warning: false
#| out-width: "70%"

dataset <- read.csv("data/glm/data_chaleur.csv")
dataset$prt_hot01 <- dataset$prt_hot/100
model <- gam(prt_hot01 ~ te(prt_veg, dist_cntr), data = dataset, family = betar)
vis.gam(model, plot.type = 'contour', type = 'response', main = '', xlab = 'Pourcentage de la surface végétalisée', ylab = 'Distance au centre-ville', too.far = 0.1, n.grid = 150)
```

Il n'y a pas de limite théorique au nombre de variables qui peuvent être ajoutées dans une *spline* d'interaction ou multivariée. Notez cependant que plus le nombre de dimensions augmente, plus la fonction à estimer est complexe et plus le volume de données nécessaire est grand et doit couvrir densément l'ensemble de l'espace d'échantillonnage multidimensionnel.

## Mise en œuvre dans R {#sec-116}

Il est possible d'ajuster des *splines* de régression dans n'importe quel *package* permettant d'ajuster des coefficients pour un modèle de régression. Il suffit de construire les bases des *splines* en amont à l'aide du *package* `splines2` et de les ajouter directement dans l'équation de régression. 
En revanche, il est nécessaire d'utiliser des *packages* spécialisés pour ajuster des *splines* de lissage. Parmi ceux-ci, `mgcv` est probablement le plus populaire du fait de sa (très) grande flexibilité, suivi des *packages* `gamlss`, `gam` et `VGAM.` Nous comparons ici les deux approches, puis nous tentons d'améliorer le modèle que nous avons ajusté pour prédire le pourcentage de surface couverte par des îlots de chaleur dans les aires de diffusion de Montréal, dans une perspective d'équité environnementale. Pour rappel, la variable dépendante est exprimée en pourcentage et nous utilisons une distribution bêta pour la modéliser.

```{r}
#| message: false
#| warning: false
library(mgcv)
# Chargement des données
dataset <- read.csv("data/gam/data_chaleur.csv", fileEncoding = "utf8")
# Ajustement du modèle de base
refmodel <- gam(hot ~
          A65Pct + A014Pct + PopFRPct + PopMVPct +
          poly(prt_veg, degree = 2)  + Arrond,
        data = dataset, family = betar(link = "logit"))
```

Dans notre première analyse de ces données, nous avons ajusté une polynomiale d'ordre 2 pour représenter un potentiel effet non linéaire de la végétation sur les îlots de chaleur. Nous remplaçons à présent ce terme par une *spline* de régression en sélectionnant quatre nœuds.

```{r}
#| message: false
#| warning: false
library(splines2)
# Création des bases de la spline
basis <- bSpline(x = dataset$prt_veg, df =4, intercept = FALSE)
# Ajouter les bases au DataFrame
basisdf <- as.data.frame(basis)
names(basisdf) <- paste('spline',1:ncol(basisdf), sep='')
dataset <- cbind(dataset, basisdf)
# Ajuster le modèle
model0 <- gam(hot ~
          A65Pct + A014Pct + PopFRPct + PopMVPct +
          spline1 + spline2 + spline3 + spline4 + Arrond,
        data = dataset, family = betar(link = "logit"))

```

Nous pouvons à présent ajuster une *spline* de lissage et laisser `mgcv` déterminer son niveau de complexité.

```{r}
#| message: false
#| warning: false
# Ajustement du modèle avec une spline simple
model1 <- gam(hot ~
          A65Pct + A014Pct + PopFRPct + PopMVPct +
          s(prt_veg)  + Arrond,
        data = dataset, family = betar(link = "logit"))
```

Notez ici que la syntaxe à employer est très simple, il suffit de spécifier `s(prt_veg)` pour indiquer à la fonction `gam` que vous souhaitez ajuster une *spline* pour la variable `prt_veg`. Nous pouvons à présent comparer l'ajustement des deux modèles en utilisant la mesure de l'AIC.

```{r}
#| message: false
#| warning: false
# Comparaison des AIC
AIC(refmodel, model0, model1)
```

Nous constatons que la valeur de l'AIC du second modèle est plus réduite, indiquant un meilleur ajustement du modèle avec une *spline* de régression. Notons cependant que la différence avec la *spline* de lissage est anecdotique (deux points de l'AIC) et que nous connaissions a priori le bon nombre de nœuds à utiliser. Pour des relations plus complexes, les *splines* de lissage ont tendance à nettement mieux performer. Voyons à présent comment représenter ces trois termes non linéaires.

```{r}
#| message: false
#| warning: false
#| label: fig-gam16
#| fig-align: center
#| fig-cap: Comparaison d'une spline et d'une polynomiale
#| out-width: "75%"

# Création d'un DataFrame de prédiction dans lequel seule
# la variable prt_veg varie.
dfpred <- data.frame(
  prt_veg = seq(min(dataset$prt_veg), max(dataset$prt_veg), 0.5),
  A65Pct = mean(dataset$A65Pct),
  A014Pct = mean(dataset$A014Pct),
  PopFRPct = mean(dataset$PopFRPct),
  PopMVPct = mean(dataset$PopMVPct),
  Arrond = "Verdun"
)

# Recréation des bases de la spline de régression
# pour les nouvelles observations
nvl_bases <- data.frame(predict(basis, newx = dfpred$prt_veg))
names(nvl_bases) <- paste('spline',1:ncol(basisdf), sep='')
dfpred <- cbind(dfpred, nvl_bases)

# Définition de la fonction inv.logit, soit l'inverse de la fonction 
# de lien du modèle pour retrouver les prédictions dans l'échelle 
# originales des données
inv.logit <- function(x){exp(x)/(1+exp(x))}

# Utilisation des deux modèles pour effectuer les prédictions
predref <- predict(refmodel, newdata = dfpred, type = 'link', se.fit = TRUE)
predmod0 <- predict(model0, newdata = dfpred, type = 'link', se.fit = TRUE)
predmod1 <- predict(model1, newdata = dfpred, type = 'link', se.fit = TRUE)

# Calcul de la valeur prédite et construction des intervalles de confiance
dfpred$polypred <- inv.logit(predref$fit)
dfpred$poly025 <- inv.logit(predref$fit - 1.96 * predref$se.fit)
dfpred$poly975 <- inv.logit(predref$fit + 1.96 * predref$se.fit)

dfpred$regsplinepred <- inv.logit(predmod0$fit)
dfpred$regspline025 <- inv.logit(predmod0$fit - 1.96 * predmod0$se.fit)
dfpred$regspline975 <- inv.logit(predmod0$fit + 1.96 * predmod0$se.fit)

dfpred$splinepred <- inv.logit(predmod1$fit)
dfpred$spline025 <- inv.logit(predmod1$fit - 1.96 * predmod1$se.fit)
dfpred$spline975 <- inv.logit(predmod1$fit + 1.96 * predmod1$se.fit)

# Créer un graphique pour afficher les résultats
ggplot(dfpred) + 
  geom_ribbon(aes(x = prt_veg, ymin = poly025, ymax = poly975),
              alpha = 0.4, color = "grey") +
  geom_ribbon(aes(x = prt_veg, ymin = spline025, ymax = spline975),
              alpha = 0.4, color = "grey") +
  geom_ribbon(aes(x = prt_veg, ymin = regspline025, ymax = regspline975),
              alpha = 0.4, color = "grey") +
  geom_line(aes(y = polypred, x = prt_veg, color = "polynomiale"),
            linewidth = 1) + 
  geom_line(aes(y = regsplinepred, x = prt_veg, color = "spline de régression"),
            linewidth = 1)+
  geom_line(aes(y = splinepred, x = prt_veg, color = "spline de lissage"),
            linewidth = 1)
```

Nous constatons que les trois termes renvoient des prédictions très similaires et qu'une légère différence n'est observable que pour les secteurs avec les plus hauts niveaux de végétation (supérieurs à 75 %).

Jusqu'ici, nous utilisons l'arrondissement dans lequel est comprise chaque aire de diffusion comme une variable nominale afin de capturer la dimension spatiale du jeu de données. Puisque nous avons abordé la notion de *splines* bivariées, il serait certainement plus efficace d'en construire une à partir des coordonnées géographiques (*x*,*y*) des centroïdes des aires de diffusion. En effet, il est plus probable que la distribution des îlots de chaleur suive un patron spatial continu sur le territoire plutôt que les délimitations arbitraires des arrondissements.

```{r}
#| message: false
#| warning: false
# Ajustement du modèle avec une spline bivariée pour l'espace
model2 <- gam(hot ~
          A65Pct + A014Pct + PopFRPct + PopMVPct +
          s(prt_veg)  + s(X,Y),
        data = dataset, family = betar(link = "logit"))
```

Notez ici que l'expression `s(X,Y)` permet de créer une *spline* bivariée à partir des coordonnées (*x, y*), soit deux colonnes présentes dans le jeu de données. Ces coordonnées sont exprimées toutes deux en mètres et n'interagissent pas ensemble au sens strict, nous devons donc ajuster une *spline* bivariée. Si vous avez besoin d'ajuster une *spline* d'interaction (notamment quand les variables sont dans des unités différentes), il est nécessaire d'utiliser une autre syntaxe `te(X,Y)` ou `t2(X,Y)` faisant appel à une structure mathématique légèrement différente, soit des *tensor product smooths*.

Puisque notre modèle intègre deux *splines*, nous devons nous assurer que nous n'avons pas de problème de concurvité, ce que nous pouvons faire avec la fonction `concurvity` du *package* `mgcv`. 

```{r}
#| message: false
#| warning: false
values <- concurvity(model2, full = FALSE)

# Worst, estimation pessimiste de la concurvité
round(values$worst,3)

# Observed, estimation optimiste de la concurvité
round(values$observed,3)

# Estimate, estimation entre deux de la concurvité
round(values$estimate,3)
```

Nous pouvons ainsi constater des niveaux de concurvité tout à fait acceptables dans notre modèle. Des valeurs supérieures à 0,8 devraient être considérées comme alarmantes, surtout si elles sont reportées pour `observed` et `estimate`.

Voyons désormais, le résumé d'un modèle GAM tel que présenté dans R.

```{r}
#| message: false
#| warning: false
summary(model2)
```

La première partie du résumé comprend les résultats pour les effets fixes et linéaires du modèle. Ils s'interprètent comme pour ceux d'un GLM classique. La seconde partie présente les résultats pour les termes non linéaires. La valeur de *p* permet de déterminer si la *spline* a ou non un effet différent de 0. Une valeur non significative indique que la *spline* ne contribue pas au modèle. Les colonnes *edf* et *Ref.df* indiquent la complexité de la *spline* et peuvent être considérées comme une approximation du nombre de nœuds. Dans notre cas, la *spline* spatiale (`s(X,Y)`) est environ 5 fois plus complexe que la *spline* ajustée pour la végétation (`s(prt_veg)`). Cela n'est pas surprenant puisque la dimension spatiale (*spline* bivariée) du phénomène est certainement plus complexe que l’effet de la végétation. Notez ici que des valeurs *edf* et *Ref.df* proches de 1 signaleraient que l'effet d’un prédicteur est essentiellement linéaire et qu'il n'est pas nécessaire de recourir à une *spline* pour cette variable.

La dernière partie du résumé comprend deux indicateurs de qualité d'ajustement, soit le R^2^ ajusté et la part de la déviance expliquée.

```{r}
#| message: false
#| warning: false
AIC(refmodel, model1, model2)
```

Nous pouvons constater que le fait d'introduire la *spline* spatiale dans le modèle contribue à réduire encore la valeur de l'AIC, et donc à améliorer le modèle. À ce stade, nous pourrions tenter de forcer la *spline* à être plus complexe en augmentant le nombre de nœuds.

```{r}
#| message: false
#| warning: false
# Augmentation de la complexité de la spline spatiale
model3 <- gam(hot ~
          A65Pct + A014Pct + PopFRPct + PopMVPct +
          s(prt_veg)  + s(X,Y,k = 40),
        data = dataset, family = betar(link = "logit"))

AIC(refmodel, model1, model2, model3)
```

Cela a pour effet d'améliorer de nouveau le modèle. Pour vérifier si l'augmentation du nombre nœuds est judicieuse, il est possible de représenter le résultat des deux *splines* précédentes. Pour ce faire, nous proposons de calculer les valeurs prédites de la *spline* pour chaque localisation dans notre terrain d'étude, en le découpant préalablement en pixels de 100 de côté. Pour cette prédiction, nous maintenons toutes les autres variables à leur moyenne respective afin d'évaluer uniquement l'effet de la *spline* spatiale.

```{r}
#| label: fig-gam17
#| message: false
#| warning: false
#| fig-align: center
#| fig-cap: Comparaison de deux splines spatiales
#| out-width: "100%"
library(viridis)
library(metR) # pour placer des étiquettes sur les isolignes

# Création d'un DataFrame fictif pour les prédictions
dfpred <- expand.grid(
  prt_veg =mean(dataset$prt_veg),
  A65Pct = mean(dataset$A65Pct),
  A014Pct = mean(dataset$A014Pct),
  PopFRPct = mean(dataset$PopFRPct),
  PopMVPct = mean(dataset$PopMVPct),
  X = seq(min(dataset$X), max(dataset$X),100),
  Y = seq(min(dataset$Y), max(dataset$Y),100)
)

dfpred$predicted1 <- predict(model2, newdata = dfpred, type = 'response')
dfpred$predicted2 <- predict(model3, newdata = dfpred, type = 'response')

# Centrage des prédictions
dfpred$predicted1 <- dfpred$predicted1 - mean(dfpred$predicted1)
dfpred$predicted2 <- dfpred$predicted2 - mean(dfpred$predicted2)

# Représentation des splines
plot1 <- ggplot(dfpred) + 
  geom_raster(aes(x = X, y = Y, fill = predicted1)) + 
  geom_point(aes(x = X, y = Y),
             size = 0.2, alpha = 0.4,
             color = "black", data = dataset)+
  geom_contour(aes(x = X, y = Y, z = predicted1), binwidth = 0.1, 
               color = "white", linetype = 'dashed') + 
  geom_text_contour(aes(x = X, y = Y, z = predicted1), 
                    color = "white", binwidth = 0.1)+
  scale_fill_viridis() +
  coord_cartesian() + 
  theme(axis.title= element_blank(),
        axis.text = element_blank(),
        axis.ticks =  element_blank()
        ) + 
  labs(subtitle = "spline de base", fill = "prédictions")

plot2 <- ggplot(dfpred) + 
  geom_raster(aes(x = X, y = Y, fill = predicted2)) + 
  geom_point(aes(x = X, y = Y),
             size = 0.2, alpha = 0.4,
             color = "black", data = dataset)+
  geom_contour(aes(x = X, y = Y, z = predicted2),
               binwidth = 0.1, color = "white", linetype = 'dashed') + 
  geom_text_contour(aes(x = X, y = Y, z = predicted2), color = "white", binwidth = 0.1)+
  scale_fill_viridis() +
  coord_cartesian()+ 
  theme(axis.title= element_blank(),
        axis.text = element_blank(),
        axis.ticks =  element_blank()
        ) + 
  labs(subtitle = "spline plus complexe", fill = "prédictions")

ggarrange(plot1, plot2, nrow = 1, ncol = 2, common.legend = TRUE, legend = 'bottom')
```

Or, il s'avère que les deux *splines* spatiales sont très similaires (@fig-gam17). Par conséquent, il est vraisemblablement plus pertinent de conserver la plus simple des deux. Notez que le Mont-Royal, compris dans le cercle central avec une isoligne à 0, est caractérisé par des valeurs plus faibles d'îlots de chaleur, alors que les quartiers centraux situés un peu plus au nord sont au contraire marqués par des pourcentages d'îlots de chaleur supérieurs de 20 points de pourcentage en moyenne.

## GAMM {#sec-117}

Bien entendu, il est possible de combiner les modèles généralisés additifs (GAM) avec les modèles à effet mixtes (GLMM) abordés dans les sections précédentes. Ces modèles généralisés additifs à effets mixtes (GAMM) peuvent facilement être mis en œuvre avec `mgcv`.

#### GAMM et interceptes aléatoires {#sec-1171}

Pour définir des constantes aléatoires, il suffit d'utiliser la notation `s(var, bs = 're')` avec `var` une variable nominale. Reprenons l'exemple précédent, mais avec cette fois-ci les arrondissements comme un intercepte aléatoire.

```{r}
#| message: false
#| warning: false
dataset$Arrond <- as.factor(dataset$Arrond)
model4 <- gam(hot ~
          A65Pct + A014Pct + PopFRPct + PopMVPct +
          s(prt_veg)  + s(Arrond, bs = "re"),
        data = dataset, family = betar(link = "logit"))
```

L'enjeu est ensuite d'extraire la variance propre à cet effet aléatoire ainsi que les valeurs des interceptes pour chaque arrondissement.

```{r}
#| message: false
#| warning: false
gam.vcomp(model4)
```

Nous constatons donc que l'écart-type de l'effet aléatoire des arrondissements est de 0,39, ce qui signifie que les effets de chaque arrondissement seront compris à 95 % entre -1,17 et 1,17 (`1.17 = 3*0.39`) sur l'échelle du prédicteur linéaire. En effet, rappelons que les effets aléatoires sont modélisés comme des distributions normales et que 95 % de la densité d'une distribution normale se situe entre -3 et +3 écarts-types. Pour extraire les interceptes spécifiques de chaque arrondissement, nous pouvons utiliser la fonction `get_random` du *package* `itsadug.`

```{r}
#| label: fig-randomconstGAM
#| fig-cap: Constantes aléatoires pour les arrondissements
#| message: false
#| warning: false
library(itsadug)
values <- get_random(model4)[[1]]
df <- data.frame(
  ri = as.numeric(values),
  arrond = names(values)
)

ggplot(df) + 
  geom_point(aes(x = ri, y = reorder(arrond, ri))) + 
  geom_vline(xintercept = 0, color = "red") + 
  labs(y = "Arrondissement", x = "intercepte aléatoire")
```

Nous constatons ainsi, à la @fig-randomconstGAM, que pour une partie des arrondissements, la densité d'îlot de chaleur est systématiquement supérieure à la moyenne régionale représentée ici par la ligne rouge (0 = effet moyen pour tous les arrondissements). Il convient alors d’améliorer ce graphique en ajoutant le niveau d'incertitude associé à chaque intercepte. Pour ce faire, nous utilisons la fonction `extract_random_effects` du *package* `mixedup.` Notez que ce *package* n'est actuellement pas disponible sur CRAN et doit être téléchargé sur **github** avec la commande suivante :

```{r}
#| message: false
#| warning: false
#| eval: false
remotes::install_github('m-clark/mixedup')
```

Avec la version 4.0.1 de R, nous avons rencontré des difficultés pour installer `mixedup`. Nous avons donc  simplement récupéré le code source de la fonction et l'avons enregistré dans un fichier de code séparé que nous appelons ici.

```{r}
#| message: false
#| warning: false
source("code_complementaire/gam_functions.R")
```

Nous pouvons ensuite procéder à l'extraction des effets aléatoires et les représenter à nouveau (@fig-randomconstGAM2).

```{r}
#| label: fig-randomconstGAM2
#| message: false
#| warning: false
#| fig-cap: Constantes aléatoires pour les arrondissements avec intervalle de confiance
df_re <- extract_random_effects.gam(model4, re = "Arrond")

ggplot(df_re) + 
  geom_errorbarh(aes(xmin = lower_2.5, xmax = upper_97.5, y = reorder(group, value))) +
  geom_point(aes(x = value, y = reorder(group, value))) + 
  geom_vline(xintercept = 0, color = "red") + 
  labs(y = "Arrondissement", x = "Intercepte aléatoire")
```

Cela permet de distinguer quels écarts sont significativement différents de 0 au seuil de 95 %. À titre de rappel, pour être significatif à ce seuil, un intervalle représenté par une ligne noire horizontale ne doit pas intersecter la ligne rouge verticale. Puisque nous utilisons ici la distribution bêta et une fonction de lien logistique, nous devons utiliser des prédictions pour simplifier l'interprétation des coefficients. Nous fixons ici toutes les variables à leur moyenne respective, sauf l'arrondissement, et calculons les prédictions dans l'échelle originale (0 à 1).

```{r}
#| label: fig-randomconstGAM3
#| message: false
#| warning: false
#| fig-cap: Prédictions pour les différents arrondissements pour une AD fictive moyenne
dfpred <- data.frame(
  A65Pct = mean(dataset$A65Pct),
  A014Pct = mean(dataset$A014Pct),
  PopFRPct = mean(dataset$PopFRPct),
  PopMVPct = mean(dataset$PopMVPct),
  prt_veg = mean(dataset$prt_veg),
  Arrond = as.character(unique(dataset$Arrond))
)

# Calculer les prédictions pour le prédicteur linéaire
dfpred$preds <- predict(model4, newdata = dfpred, type = "link")

# Calculer l'intervalle de confiance en utilisant les valeurs
# extraites avec extract_random_effects
dfpred <- dfpred[order(dfpred$Arrond),]
df_re <- df_re[order(df_re$group),]

dfpred$lower <- dfpred$preds - 1.96*df_re$se
dfpred$upper <- dfpred$preds + 1.96*df_re$se

# Il nous reste juste à reconvertir le tout dans l'unité d'origine
# en utilisant l'inverse de la fonction logistique
inv.logit <- function(x){exp(x)/(1+exp(x))}

dfpred$lower <- inv.logit(dfpred$lower)
dfpred$upper <- inv.logit(dfpred$upper)
dfpred$preds <- inv.logit(dfpred$preds)

ggplot(dfpred) + 
  geom_errorbarh(aes(xmin = lower, xmax = upper, y = reorder(Arrond, preds))) +
  geom_point(aes(x = preds, y = reorder(Arrond, preds))) + 
  geom_vline(xintercept = mean(dfpred$preds), color = "red") + 
  labs(y = "Arrondissement", x = "intercepte aléatoire")

```

Nous constatons ainsi, à la @fig-randomconstGAM3, que pour une hypothétique aire de diffusion moyenne, la différence de densité d'îlot de chaleur peut être de 0,32 (32 % de la surface de l'AD) entre les arrondissements Verdun et Dollard-des-Ormeaux.

#### GAMM et coefficients aléatoires {#sec-1172}

En plus des interceptes aléatoires, il est aussi possible de définir des coefficients aléatoires. Reprenons notre exemple et tentons de faire varier l'effet de la variable `PopFRPct` en fonction de l'arrondissement.

```{r}
#| message: false
#| warning: false
model5 <- gam(hot ~
          A65Pct + A014Pct + PopFRPct + PopMVPct +
          s(prt_veg)  + s(Arrond, bs = "re") + 
          s(PopFRPct, Arrond, bs = "re"),
        data = dataset, family = betar(link = "logit"))
```

Notez ici une distinction importante! Le modèle n'assume aucune corrélation entre les coefficients aléatoires pour la variable `PopFRPct` et pour les constantes aléatoires. Il est présumé que ces deux effets proviennent de deux distributions normales distinctes. En d'autres termes, le modèle ne dispose pas des paramètres nécessaires pour vérifier si les arrondissements avec les constantes les plus fortes (avec des densités supérieures d'îlot de chaleur) sont aussi des arrondissements dans lesquels l'effet de la variable `PopFRPct` est plus prononcé (et vice-versa). Pour plus d'informations sur cette distinction, référez-vous à la [section @sec-0924].

```{r}
#| message: false
#| warning: false
AIC(model4, model5)
```

Ce dernier modèle présente une valeur de l'AIC plus faible et serait donc ainsi mieux ajusté que notre modèle avec seulement un intercepte aléatoire. Nous pouvons donc extraire les coefficients aléatoires et les représenter à la @fig-randomconstcoefGAM.

```{r}
#| label: fig-randomconstcoefGAM
#| message: false
#| warning: false
#| fig-cap: Pentes et constantes aléatoires pour les arrondissements
df_re <- extract_random_effects.gam(model5)
df_re <- subset(df_re, df_re$effect == 'PopFRPct')

ggplot(df_re) + 
  geom_errorbarh(aes(xmin = lower_2.5, xmax = upper_97.5, y = reorder(group, value))) +
  geom_point(aes(x = value, y = reorder(group, value))) + 
  geom_vline(xintercept = 0, color = "red") + 
  labs(y = "Arrondissement", x = "coefficient aléatoire")
```

Nous constatons notamment que seuls trois arrondissements ont des coefficients aléatoires significativement différents de 0. Ainsi, pour les arrondissements Anjou et Plateau-Mont-Royal, les coefficients aléatoires sont respectivement de `r tofr(df_re[df_re$group == "Anjou",]$value[[1]])` et `r tofr(df_re[df_re$group == "Le Plateau-Mont-Royal",]$value[[1]])`, et viennent donc se retrancher à la valeur moyenne régionale de `r tofr(round(model5$coefficients[[4]],4))` qui atteint alors presque 0. Du point de vue de l'interprétation, nous pouvons en conclure que le groupe des personnes à faible revenu ne subit pas de surexposition aux îlots de chaleur à l'échelle des AD dans ces arrondissements.

En revanche, dans l'arrondissement Mercier-Hochelaga-Maisonneuve, la situation est à l'inverse plus systématiquement en défaveur des populations à faible revenu, avec une taille d'effet près de deux fois supérieure à la moyenne régionale. En effet, l'effet moyen régional (coefficient fixe) est de `r tofr(round(model5$coefficients[[4]],4))`, auquel vient s'ajouter l'effet spécifique (coefficient aléatoire) de Mercier-Hochelaga-Maisonneuve, soit `r tofr(df_re[df_re$group == "Mercier-Hochelaga-Maisonneuve",]$value[[1]])`, pour un effet total de `r tofr(df_re[df_re$group == "Mercier-Hochelaga-Maisonneuve",]$value[[1]] + round(model5$coefficients[[4]],4))`

::: bloc_aller_loin
::: bloc_aller_loin-header
::: bloc_aller_loin-icon
:::
**Des effets aléatoires plus complexes dans les GAMM**
:::
::: bloc_aller_loin-body
Il est possible de spécifier des GAMM avec des effets aléatoires plus complexes autorisant, par exemple, des corrélations entre les différents effets / niveaux. Il faut pour cela utiliser la fonction `gamm` de `mgcv` ou la fonction `gamm4` du *package* `gamm4.` La première offre plus de flexibilité, mais la seconde est plus facile à utiliser et doit être privilégiée quand un modèle comporte un très grand nombre de groupes dans un effet aléatoire, ou lorsque la distribution du modèle n'est pas gaussienne. La fonction `gamm` permet d'ajuster des modèles non gaussiens, mais elle utilise une approche appelée PQL (*Penalized Quasi-Likelihood* en anglais) connue pour être moins stable et moins précise.

Cependant, dans l'exemple de cette section, nous utilisons un modèle GAMM avec une distribution bêta, ce qui n'est actuellement pas supporté par les fonctions `gamm` et `gamm4`. Pour un modèle GAMM plus complexe utilisant une distribution bêta, il est nécessaire d'utiliser le *package* `gamlss`, mais ce dernier utilise aussi une approche de type PQL. Nous montrons tout de même ici comment ajouter un modèle qui inclut une corrélation entre les deux effets aléatoires de l'exemple précédent. Notez ici que le terme `re` apparaissant dans la formule permet de spécifier un effet aléatoire en utilisant la syntaxe du *package* `nlme.` Plus spécifiquement, `gamlss` fait un pont avec `nlme` et utilise son algorithme d'ajustement au sein de ces propres routines. De même, le terme `pb` permet de spécifier une *spline* de lissage dans le même esprit que `mgcv.` Il est également possible d'utiliser le terme `ga` faisant le lien avec `mgcv` et de profiter de sa flexibilité dans `gamlss`.

```{r}
#| message: false
#| warning: false
#| results: hide
library(gamlss)
library(gamlss.add)

model6 <- gamlss(hot ~
          pb(prt_veg) + 
          re(fixed = ~ A65Pct + A014Pct + PopFRPct + PopMVPct, 
             random = ~(1 + PopFRPct)|Arrond),
        data = dataset, family = BE(mu.link = "logit"))
```

Nous pouvons ensuite accéder à la partie du modèle qui nous intéresse, soit celle concernant les effets aléatoires.

```{r}
#| message: false
#| warning: false
randomPart <- model6$mu.coefSmo[[2]]
print(randomPart)
```


À lecture de la partie du résumé consacrée aux résultats pour les effets aléatoires, nous constatons que la corrélation entre les interceptes aléatoires et les coefficients aléatoires est de -0,65. Cela signifie que pour les arrondissements avec des interceptes élevés (plus grande proportion d'îlots de chaleur), l'effet de la variable `PopFRPct` tend à être plus faible. Autrement dit, dans les arrondissements avec beaucoup d'îlots de chaleur, les personnes à faible revenu ont tendance à être moins exposées, tel qu'illustré à la @fig-corrrandom).

```{r}
#| label: fig-corrrandom
#| fig-cap: Relation entre les effets aléatoires des arrondissements et la variable population à faible revenu
#| fig-align: center
#| message: false
#| warning: false
#| out-width: "75%"

df <- ranef(randomPart)
df$arrond <- rownames(df)
names(df) <- c('Intercept', 'PopFRPct', 'Arrondissement')

ggplot(df) + 
  geom_hline(yintercept = 0, color = "red") +
  geom_vline(xintercept = 0, color = "red") +
  geom_point(aes(x = Intercept, y = PopFRPct))
```
:::
:::

## Quiz de révision du chapitre {#sec-118}

```{r}
#| label: quizGAM
#| echo: false
#| warning: false
#| results: asis
source("code_complementaire/QuizzFunctions.R")
quizz_gam <- quizz("quiz/Chapitre11-GAM.yml", "quizz_gam")
render_quizz(quizz_gam)
```