Librairies et jeux de données

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

• ISIT.csv

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

• ggplot2
• itsadug
• mgcv
  

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:

1. La validation d'un GAM
2. Choisir une autre distribution
3. Changer la fonction de base
4. 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.

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, 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$$

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).

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.

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:

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.

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?

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

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

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

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).

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()

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
...

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$.

GAM

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

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.

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.

Défi 1

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:

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?

Défi 1 - Solution

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

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

Défi 1 - Solution

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.

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

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

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

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

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:

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!

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.

Y a-t-il une meilleure façon de sélectionner les emplacements des nœuds?

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é.

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:

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

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).

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:

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:

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

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?

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

Degrés de liberté effectifs (EDF)

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é.

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:

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):

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):

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

GAM à plusieurs variables linéaires et lisses

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.

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) :

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) :

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.

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

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:

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

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

Défi 2

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.

Défi 2 - Solution

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.

# 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)

# 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

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.

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

Défi 2 - Solution

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

Défi 2 - Solution

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?

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

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():

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

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.

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.

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

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

Interaction: variables non-linéaire et qualitatif

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

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().

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

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

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.

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.

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

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

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?

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.

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é)

Ces situations se produisent fréquemment !

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 :

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":

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?

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, nous 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é:

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

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?

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:

smooth_interact <- smooth_interact_k60

Validation GAM: Choisir une autre distribution

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

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

Validation GAM: Choisir une autre distribution

Hétéroscédasticité marquée et points extrêmes dans les résidus

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 :

?family.mgcv

Défi 3

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 ?

Indice:

# 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

1. Premièrement, construisons un nouveau modèle avec une distribution Tweedie un lien log.

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

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:

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

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

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

Défi 3 - Solution

3. Comparer smooth_interact_tw avec smooth_interact. Lequel est meilleur?

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

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!

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() 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() 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() 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 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 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 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:

data(nottem)

# 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)

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

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.

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

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

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.

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.

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

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).

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.

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é?

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 variable qualitative utilisée pour l'effet aléatoire; x0 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.

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

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

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

# 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

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 :

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:

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

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

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.

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

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 :

plot(gamm_int_slope, select = 3)

GAMM avec une surface lisse aléatoire

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

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

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

plot(gamm_smooth, select=1)

GAMM avec une surface lisse aléatoire

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

Comparaison de GAMM

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

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 par Simon Wood (auteur de la librairie mgcv).
• Le site web de Simon Wood, maths.ed.ac.uk/~swood34/.
• Le blogue de Gavin Simpson, From the bottom of the heap.
• La librairie gratia de Gavin Simpson for GAM visualisation in ggplot2.
• Le cours Generalized Additive Models: An Introduction with R de Noam Ross.
• Overview GAMM analysis of time series data tutoriel de Jacolien van Rij.
• Hierarchical generalized additive models in ecology: an introduction with mgcv de Pedersen et al. (2019).

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

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.

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

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).

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

Qu'est ce que représente l'intercepte dans ce modèle?

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

Qu'est ce que le terme de lissage indique?

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

Que représente l'intercepte dans ce modèle?

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 :

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

Qu'est ce que le terme de lissage indique?

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$

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?

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() :

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$.