AVIS IMPORTANT

Depuis l'automne 2021, ce wiki a été discontinué et n'est plus activement développé.

Tout le matériel mis à jour et les annonces pour la série d'ateliers R du CSBQ se trouvent maintenant sur le site web de la série d'ateliers R du CSBQ. Veuillez mettre à jour vos signets en conséquence afin d'éviter les documents périmés et/ou les liens brisés.

Merci de votre compréhension,

Vos coordonnateurs de la série d’ateliers R du CSBQ.

Ateliers R du CSBQ

Cette série de 10 ateliers guide les participants à travers les étapes requises afin de maîtriser le logiciel R pour une grande variété d’analyses statistiques pertinentes en recherche en biologie et en écologie. Ces ateliers en libre accès ont été créés par des membres du CSBQ à la fois pour les membres du CSBQ et pour la grande communauté d’utilisateurs de R.

Le contenu de cet atelier a été révisé par plusieurs membres du CSBQ. Si vous souhaitez y apporter des modifications, veuillez SVP contacter les coordonnateurs actuels de la série, listés sur la page d'accueil

AVIS IMPORTANT: MISES À JOUR MAJEURES

Mise à jour de mars 2021: Ce wiki n'est plus activement développé ou mis à jour. Le matériel mis à jour pour la série d'ateliers R du CSBQ est maintenant hébergé sur la page GitHub des ateliers R du CSBQ.

Le matériel disponible inclut;

  1. La présentation Rmarkdown pour cet atelier;
  2. Un script R qui suit la présentation - *en construction*.
  3. Le matériel écrit qui accompagne la présentation en format bookdown - *en construction*.

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

Développé par : Eric Pedersen et Zofia Taranu

Révision par : Cédric Frenette Dussault

Résumé: L'objectif de l'atelier d'aujourd'hui sera d'examiner ce que nous entendons par un modèle non-linéaire et comment les GAMs (modèles additifs généralisés) nous permettent de modéliser les relations non-linéaires. Nous examinerons également comment tracer et interpréter ces relations non-linéaires, comment ajouter des interactions, comment prendre en compte la non-indépendance des données (e.g. erreurs autocorrélées) et comment inclure des effets aléatoires en se basant sur les ateliers précédents. Enfin, nous allons brièvement aborder la mécanique derrière le fonctionnement des GAMs.

Lien vers la nouvelle présentation Rmarkdown

Lien vers la présentation Prezi associée : Prezi

Télécharger le script R et les données pour cette leçon:

  1. Script
  2. Données

Objectifs de l'atelier

  1. Utiliser la librairie mgcv pour modéliser les relations non linéaires
  2. Évaluer la sortie d'un 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 la raison pour laquelle ç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

Prérequis: Expérience du logiciel R (assez pour être en mesure d'exécuter un script et d'examiner les données et les objets dans R) et une connaissance de base de la régression simple (vous devez savoir ce qu'on entend par une régression linéaire et une ANOVA).

0. Le modèle linéaire ... et où il échoue

Qu'est-ce qu'un modèle linéaire ? La régression linéaire est ce que la plupart des gens apprennent avant tout en statistiques et est parmi les méthodes les plus performantes. Elle nous permet de modéliser une variable réponse en fonction de facteurs prédictifs et d'une erreur résiduelle. Tel que vu dans l'atelier sur les modèles linéaires, la régression fait cependant quatre suppositions importantes :

  1. l'erreur est distribuée normalement
  2. la variance des erreurs est constante
  3. chaque erreur est indépendante des autres (homoscédasticité)
  4. la réponse est linéaire: {y} = {β_0} + {β_1}{x}

Il n'y a qu'une façon pour qu'un modèle linéaire soit correctement appliqué :

et pourtant tant de façons pour qu'il ne le soit pas :

Alors, comment résoudre ce problème ? Nous devons premièrement savoir ce que le modèle de régression cherche à faire :

  • ajuster une ligne qui passe au milieu des données,
  • faire cela sans sur-ajuster les données, c'est-à-dire en passant une ligne entre chaque point.

Les modèles linéaires le font en trouvant la meilleure ligne droite qui passe à travers les données. En revanche, les modèles additifs font cela en ajustant une courbe à travers les données, mais tout en contrôlant le degré de courbure de la ligne (plus d'information sur cela plus bas).

1. Introduction aux GAMs

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

| Générer et tracer des données
library(ggplot2)
set.seed(10)
n = 250
x = runif(n,0,5)
y_model = 3*x/(1+2*x)
y_obs = rnorm(n,y_model,0.1)
data_plot = qplot(x, y_obs) + 
            geom_line(aes(y=y_model)) + 
            theme_bw()
print(data_plot)

Si nous modélisions cette relation par une régression linéaire, les résultats ne respecteraient pas les suppositions énumérées ci-dessus. Commençons par modéliser une régression en utilisant la méthode des moindres carrés en utilisant la fonction gam() de la librairie mgcv - donc en tant que modèle linéaire (nous verrons plus bas comment utiliser cette fonction pour spécifier un terme non linéaire).

| Modèle linéaire
library(mgcv)
linear_model = gam(y_obs~x)
model_summary=summary(linear_model)
print(model_summary)
data_plot = data_plot+
             geom_line(colour="red",
             aes(y=fitted(linear_model)))
print(data_plot)

Nous pouvons constater à partir du sommaire que notre modèle linéaire explique une grande partie de la variance (R2(adj) = 0.639). Toutefois, les graphiques de diagnostic des résidus du modèle montrent que l'écart type ne suit pas une distribution normale et que la variance n'est pas homoscédastique. De plus, il reste un patron non-linéaire important. Essayons maintenant de résoudre ce problème en ajustant les données avec un terme non linéaire.

Nous reviendrons sur ceci un peu plus tard, mais brièvement, les GAMs sont une forme non paramétrique de la régression où le beta*xi d'une régression linéaire est remplacé par une fonction de lissage des variables explicatives, f (xi), et le modèle devient :

{y_i} = f(x_{i}) + {ε_i}

où yi est la variable réponse, xi est la covariable, et f est la fonction lissage.

Étant donné que la fonction de lissage f(xi) est non linéaire et locale, l'ampleur de l'effet de la variable explicative peut varier en fonction de la relation entre la variable et la réponse. Autrement dit, contrairement à un coefficient fixe beta, la fonction f peut changer tout au long du gradient xi. Le degré de lissage de f est contrôlée en utilisant une régression pénalisée qui est déterminée automatiquement à l'aide d'une méthode de validation croisée généralisée (GCV) de la librairie mgcv (Wood 2006).

Avec gam() les termes non linéaires sont spécifiés par des expressions de la forme: s(x).

| GAM
gam_model = gam(y_obs~s(x))
summary(gam_model)
data_plot = data_plot +  
     geom_line(colour="blue",aes(y=fitted(gam_model)))
print(data_plot)

La variance expliquée par notre modèle a augmenté de 20% (R2(adj) = 0.859) et quand on compare l'ajustement du modèle linéaire (rouge) au modèle non-linéaire (bleu), il est clair que l'ajustement de ce dernier est relativement meilleur.

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

| smooth plot
plot(gam_model)

Comment utilisons-nous les GAMs pour savoir si un modèle linéaire est suffisant pour modéliser nos données ? Nous pouvons utiliser les fonctions gam() et anova() pour tester formellement si une hypothèse de linéarité est justifiée. Nous devons simplement le configurer de sorte que notre modèle non-linéaire est emboîté dans notre modèle linéaire; c'est à dire, nous devons créer un objet qui inclut à la fois x (linéaire) et s(x) (non-linéaire). En utilisant la fonction anova(), on vérifie si l'ajout de s(x) au modèle avec seulement x comme covariable est justifié par les données.

| Test de linéarité
linear_model = gam(y_obs~x)
nested_gam_model = gam(y_obs~s(x)+x)
print(anova(linear_model, nested_gam_model, test="Chisq"))

Le terme non linéaire est significatif: 

 Analysis of Deviance Table
 Model 1: y_obs ~ x
 Model 2: y_obs ~ s(x) + x
       Resid. Df    Resid. Dev     Df        Deviance    Pr(>Chi)    
 1    248.00        6.5846                              
 2    240.68        2.4988         7.3168    4.0858      < 2.2e-16 ***
 ---
 Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

DÉFI 1

Nous allons maintenant essayer cela avec d'autres données générées aléatoirement. Nous allons d'abord générer les données. Ensuite, nous allons ajuster un modèle linéaire et un GAM à la relation entre x_test et y_test_obs. Quels sont les degrés de libertés effectifs du terme non-linéaire ? Déterminez si l'hypothèse de linéarité est justifiée pour ces données.

Générer des données
n <- 250
x_test <- runif(n,-5,5)
y_test_fit <- 4*dnorm(x_test)
y_test_obs <- rnorm(n,y_test_fit, 0.2)

Réponse au défi 1

2. Plusieurs termes 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 un ensemble de données générées automatiquement par mgcv.

| Générer des données non-linéaires
gam_data = gamSim(eg=5)
head(gam_data)

Nous allons voir comment nous pouvons prédire la variable réponse y en fonction des autres variables. Commençons par un modèle de base comprenant un terme non linéaire (x1) et un facteur qualitatif (X0 avec 4 niveaux).

| Générer des données non linéaires 
basic_model = gam(y~x0+s(x1), data= gam_data)
basic_summary = summary(basic_model)
print(basic_summary$p.table)
print(basic_summary$s.table)
plot(basic_model)

Ici, la sortie de p.table fournit le tableau de résultats pour chaque terme paramétrique et le tableau s.table nous donne les résultats du terme non linéaire. Notez que pour le second tableau, la courbure du terme non linéaire s(X1) est indiquée par le paramètre edf (degrés de libertés effectifs); plus la valeur de l'edf est élevée, plus la non-linéarité est forte. Une valeur élevée (8-10 ou plus) signifie que la courbe est fortement non linéaire, alors qu'une courbe avec un edf égal à 1 est une ligne droite. En revanche, dans la régression linéaire, les degrés de libertés du modèle sont équivalents au nombre de paramètres libres non redondants p dans le modèle (et les degrés de libertés résiduels sont égaux à n-p). Nous reviendrons plus tard sur le concept d'edf. 

 print(basic_summary$p.table)
             Estimate Std. Error   t value     Pr(>|t|)
 (Intercept) 8.550030  0.3655849 23.387258 1.717989e-76
 x02         2.418682  0.5165515  4.682364 3.908046e-06
 x03         4.486193  0.5156501  8.700072 9.124666e-17
 x04         6.528518  0.5204234 12.544629 1.322632e-30
 > print(basic_summary$s.table)
            edf   Ref.df        F      p-value
 s(x1) 1.923913 2.406719 42.84268 1.076396e-19

Dans notre modèle de base, l'edf du terme non linéaire s(x1) est ~ 2, ce qui indique une courbe non linéaire. Le graphique du modèle illustre bien la forme de ce terme non linéaire :

Nous pouvons ajouter un second terme x2, mais spécifier une relation linéaire avec Y (i.e. les GAMs peuvent inclure à la fois des termes linéaires et non linéaires dans le même modèle). Ce nouveau terme linéaire x2 sera présenté dans le tableau p.table, pour lequel une estimation du coefficient de régression sera indiquée. Dans le tableau s.table, nous retrouvons encore une fois le terme non linéaire s(x1) et son paramètre de courbure.

| deux termes
two_term_model <- gam(y~x0+s(x1)+x2, data=gam_data)
two_term_summary <- summary(two_term_model)
print(two_term_summary$p.table)
print(two_term_summary$s.table)

Pour évaluer si la relation entre y et x2 est non linéaire, on peut modéliser x2 avec une fonction non linéaire. Tel que vu auparavant, nous pouvons utiliser une ANOVA pour tester si le terme non linéaire est nécessaire.

| Deux termes non linéaires 
two_smooth_model <- gam(y~x0+s(x1)+s(x2), data=gam_data)
two_smooth_summary <- summary(two_smooth_model)
print(two_smooth_summary$p.table)
print(two_smooth_summary$s.table)
plot(two_smooth_model,page=1)

Lorsqu'il y a plus d'une variable d'incluse dans le modèle, comme ci-dessus, la réponse ajustée peut-être partitionnée entre les contributions de chaque variable. Ici, nous pouvons évaluer l'effet de chaque variable où l'axe des ordonnées représente la contribution (effet) de chaque covariable à la réponse ajustée centrée sur 0. Si l'intervalle de confiance chevauche zéro pour certaines valeurs de x, cela indique que l'effet est non significatif. Lorsque la contribution varie selon l'axe x, un changement de cette variable cause un changement de la variable réponse.

| ANOVA
anova(basic_model,two_term_model,two_smooth_model, test="Chisq")
  Analysis of Deviance Table
  Model 1: y ~ x0 + s(x1)
  Model 2: y ~ x0 + s(x1) + x2
  Model 3: y ~ x0 + s(x1) + s(x2)
    Resid. Df Resid. Dev      Df Deviance  Pr(>Chi)    
  1    394.08     5231.6                               
  2    393.10     4051.3 0.97695   1180.2 < 2.2e-16 ***
  3    385.73     1839.5 7.37288   2211.8 < 2.2e-16 ***
  ---
  Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Le meilleur modèle est le modèle avec deux fonctions non linéaires pour x1 et pour x2.

DÉFI 2

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

Réponse au défi 2

3. Interactions

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

  • si une variable est quantitative et l'autre est qualitative, on utilise l'argument by → s(x, by=facteur),
  • si les deux variables sont quantitatives, on inclut les deux termes sous une même fonction non linéaire → s(x1, x2).

L'argument by permet de faire varier un terme non linéaire selon les différents niveaux d'un facteur. Nous allons examiner ceci en utilisant notre variable qualitative x0 et examiner si la non-linérité de s(x2) varie selon les différents niveaux de x0. Pour déterminer si les courbes diffèrent significativement entre les niveaux du facteur, nous allons utiliser une ANOVA sur l'interaction.

| Interaction qualitative 
categorical_interact <- gam(y~x0+s(x1)+s(x2,by=x0),data=gam_data)
categorical_interact_summary <- summary(categorical_interact)
print(categorical_interact_summary$s.table)
plot(categorical_interact,page=1)
# ou nous pouvons utiliser la fonction vis.gam où theta représente la rotation du plan x-y
vis.gam(categorical_interact,view=c("x2","x0"),theta=40,n.grid=500,border=NA) 
anova(two_smooth_model, categorical_interact,test="Chisq")

Nous pouvons constater à partir du graphique que les formes des termes non linéaires sont comparables entre les quatre niveaux de x0. L'ANOVA le confirme également (déviance = 98,6, p = 0,2347).

Ensuite, nous allons examiner l'interaction non linéaire entre deux termes quantitatifs, x1 et x2. Cette fois-ci, l'argument by est supprimé.

| Interaction quantitative
smooth_interact <- gam(y~x0+s(x1,x2),data=gam_data)
smooth_interact_summary <- summary(smooth_interact)
print(smooth_interact_summary$s.table)
plot(smooth_interact,page=1,scheme=3)
# plot(smooth_interact,page=1,scheme=1) donne un graphique comparable à vis.gam()
vis.gam(smooth_interact,view=c("x1","x2"),theta=40,n.grid=500,border=NA) 
anova(two_smooth_model,smooth_interact,test="Chisq")

L'interaction entre s(x1) et s(x2) est significative et le graphique en deux dimensions illustre très bien cette interaction non linéaire. La relation entre y et x1 change en fonction de la valeur de x2. Vous pouvez changez la valeur de l'argument theta pour tourner l'axe du graphique. Si vous prévoyez exécuter un grand nombre de graphiques, supprimez l'argument n.grid = 500, car ceci fait appel à des calculs intensifs et ralentit R.

4. Changer la fonction de base

Sans entrer dans le détail, sachez qu'il est possible de modifier le modèle de base que nous avons vu avec :

  • des fonctions plus complexes en modifiant la fonction de base (par exemple, cyclique),
  • d'autres distributions : tout ce que vous pouvez faire avec un GLM (tel que spécifier l'argument family) est possible avec les GAMs,
  • des modèles à effets mixtes en utilisant la fonction gamm ou la fonction gamm4 de la librairie gamm4.

Nous allons d'abord jeter un coup d’œil au changement de la fonction de base puis une introduction rapide aux autres distributions et les GAMMs (modèles additifs généralisés à effets mixtes) suivra. Commençons par regarder un cas où modifier la fonction de base peut être utile, soit avec des données cycliques.

Imaginez que vous avez une série temporelle de données climatiques, divisées en mesures mensuelles, et que vous voulez déterminer s'il y a une tendance de température annuelle. Nous allons utiliser la série temporelle de température de Nottingham pour cette section :

| Fonction de base cyclique
data(nottem)
n_years <- length(nottem)/12
nottem_month <- rep(1:12, times=n_years)
nottem_year <- rep(1920:(1920+n_years-1),each=12)
nottem_plot <- qplot(nottem_month,nottem, 
                    colour=factor(nottem_year), 
                    geom="line") + theme_bw()
print(nottem_plot)

En utilisant les données nottem, nous avons créé trois nouveaux vecteurs ; n_years correspond au nombre d'années de données (20 ans), nottem_month est un 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) et nottem_year est une variable où l'année correspondant à chaque mois est fournie.

Données mensuelles des années 1920 à 1940:

Pour modéliser cela, nous devons utiliser ce qu'on appelle un “spline cubique cyclique”, ou cc, pour modéliser les effets du mois et de l'année.

| gam cyclique 
year_gam <- gam(nottem~s(nottem_year)+s(nottem_month, bs="cc"))
summary(year_gam)$s.table
plot(year_gam,page=1, scale=0)

Il y a une hausse d'environ 1 - 1,5ºC au cours de la série, mais au cours d'une année, il y a une variation d'environ 20ºC. Les données réelles varient autour de ces valeurs prédites et ceci représente donc la variance inexpliquée. Ici, nous pouvons voir l'un des avantages très intéressants de l'utilisation des GAMs. Nous pouvons soit tracer la surface réponse (valeurs prédites) ou les termes (contribution de chaque covariable) tel qu'indiqué ci-haut. Vous pouvez imaginer ce dernier en tant qu'une illustration de la variation des coefficients de régression et comment leur contribution (ou taille de leur effet) varie au fil du temps. Dans le premier graphique, nous voyons que les contributions positives de la température sont survenues après 1930.

Sur des échelles de temps plus longues, en utilisant par exemple des données paléolimnologiques, d'autres (Simpson & Anderson 2009;. Fig 3c) ont utilisé des GAMs pour tracer la contribution (effet) de la température sur la composition d'algues dans les lacs afin d'illustrer comment les contributions significatives ont seulement eu lieu au cours de deux périodes extrêmement froides (c'est-à-dire, la contribution est importante lorsque les intervalles de confiance ne recoupent pas la valeur de zéro à environ 300 et 100 ans AVJC). Cela a permis aux auteurs de non seulement déterminer combien de variance est expliquée par la température au cours des derniers siècles, mais aussi de repérer dans le temps cet effet significatif. Si cela vous intéresse, le code pour tracer soit la surface de réponse (type = “response”) ou les termes (type = “terms”) est disponible ci-dessous. Lorsque les termes sont sélectionnés, vous obtiendrez la même figure que celle ci-dessus.

Graphique de contribution vs réponse ajustée

5. Intro rapide aux GAMMs

La non-indépendance des données

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

  • une structure de corrélation pour modéliser les résidus autocorrélés (autorégressif (AR), moyenne mobile (MA), ou une combinaison des deux (ARMA))
  • des effets aléatoires qui modélisent l'indépendance entre les observations d'un même site.

En plus de changer la fonction de base, nous pouvons aussi complexifier le modèle en intégrant une structure d'auto-corrélation (ou même des effets mixtes) en utilisant les fonctions gamm ou gamm4. Pour commencer, nous allons jeter un coup d’œil au premier cas ; un modèle avec autocorrélation temporelle dans les résidus. Ré-examinons le modèle de la température de Nottingham; nous allons vérifier si les résidus sont corrélés en faisant appel à la fonction (partielle) d'autocorrélation.

Erreurs corrélées
par(mfrow=c(1,2))
acf(resid(year_gam), lag.max = 36, main = "ACF")
pacf(resid(year_gam), lag.max = 36, main = "pACF")

Les graphiques des fonctions d'autocorrélation suggèrent qu'un modèle AR de faible ordre est nécessaire (avec un ou deux intervalles de temps décalés), donc nous pouvons évaluer deux modèles; ajouter un AR(1) ou un AR(2) au modèle de la température de Nottingham et évaluer le meilleur avec une ANOVA.

modèles AR
year_gam <- gamm(nottem~s(nottem_year)+s(nottem_month, bs="cc"))
year_gam_AR1 <- gamm(nottem~s(nottem_year)+s(nottem_month, bs="cc"),
                     correlation = corARMA(form = ~ 1|nottem_year, p = 1))
year_gam_AR2 <- gamm(nottem~s(nottem_year)+s(nottem_month, bs="cc"),
                     correlation = corARMA(form = ~ 1|nottem_year, p = 2))
anova(year_gam$lme,year_gam_AR1$lme,year_gam_AR2$lme)
                 Model df      AIC      BIC    logLik   Test   L.Ratio p-value
year_gam$lme         1  5 1109.908 1127.311 -549.9538                         
year_gam_AR1$lme     2  6 1101.218 1122.102 -544.6092 1 vs 2 10.689206  0.0011
year_gam_AR2$lme     3  7 1101.598 1125.962 -543.7988 2 vs 3  1.620821  0.2030

Le modèle avec la structure AR(1) prévoit une augmentation significative comparativement au premier modèle (LRT = 10,69, p = 0,0011), mais il y a très peu d'intérêt à considérer le modèle AR(2) (LRT = 1,62, b = 0,203).

Modélisation avec effets mixtes

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 représente une variable qualitative utilisée pou l'effet aléatoire et x0 est un 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.

Nous examinerons d'abord un GAMM avec un interception aléatoire. Tel que vu précédemment, nous allons utiliser gamSim() pour générer un ensemble de données, cette fois-ci avec une composante d'effet aléatoire. Ensuite, nous construirons un modèle avec un intercepte aléatoire en utilisant fac comme facteur aléatoire.

Intercepte aléatoire 
# Générez des données
gam_data2 <- gamSim(eg=6)
head(gam_data2)
 
# Faites rouler un modèle avec intercepte aléatoire
gamm_intercept <- gam(y ~ s(x0) + s(fac, bs="re"), data=gam_data2)
summary(gamm_intercept)

Notez le terme aléatoire dans le tableau. Vous pouvez le visualiser:

graphique intercepte aléatoire
plot(gamm_intercept, select=2) 
# select=2 parce que le terme aléatoire se trouve sur la 2e ligne du tableau sommaire.

Une fonction de traçage vraiment intéressante que nous allons maintenant utiliser est le plot_smooth de la librairie itsadug. Contrairement au graphique par défaut plot.gam, cette fonction présente l'effet additionné du GAMM avec l'option de ne pas inclure les courbes aléatoires dans le graphique. Ici, nous allons premièrement tracer l'effet combiné de x0 (sans les niveaux de l'effet aléatoire) et ensuite une courbe pour les quatre niveaux de fac:

| Graphique des 4 niveaux du facteur fac
par(mfrow=c(1,2), cex=1.1)
plot_smooth(gamm_intercept, view="x0", rm.ranef=TRUE, main="intercept + s(x1)", rug=FALSE)
plot_smooth(gamm_intercept, view="x0", cond=list(fac="1"), 
            main="... + s(fac)", col='orange', ylim=c(8,21), rug=FALSE)
plot_smooth(gamm_intercept, view="x0", cond=list(fac="2"), add=TRUE, col='red')
plot_smooth(gamm_intercept, view="x0", cond=list(fac="3"), add=TRUE, col='purple')
plot_smooth(gamm_intercept, view="x0", cond=list(fac="4"), add=TRUE, col='turquoise')

Ensuite, nous allons générer et tracer un modèle avec une pente aléatoire :

| pente aléatoire 
gamm_slope <- gam(y ~ s(x0) + s(x0, fac, bs="re"), data=gam_data2)
summary(gamm_slope)
 
plot_smooth(gamm_slope, view="x0", rm.ranef=TRUE, main="intercept + s(x0)", rug=FALSE)
plot_smooth(gamm_slope, view="x0", cond=list(fac="1"), 
            main="... + s(fac)", col='orange',ylim=c(7,22), rug=FALSE)
plot_smooth(gamm_slope, view="x0", cond=list(fac="2"), add=TRUE, col='red')
plot_smooth(gamm_slope, view="x0", cond=list(fac="3"), add=TRUE, col='purple')
plot_smooth(gamm_slope, view="x0", cond=list(fac="4"), add=TRUE, col='turquoise')

Nous allons maintenant inclure à la fois un intercepte et une pente aléatoires.

Intercepte et pente aléatoires 
gamm_int_slope <- gam(y ~ s(x0) + s(fac, bs="re") 
                      + s(fac, x0, bs="re"), data=gam_data2)
summary(gamm_int_slope)
 
plot_smooth(gamm_int_slope, view="x0", rm.ranef=TRUE, main="intercept + s(x0)", rug=FALSE)
plot_smooth(gamm_int_slope, view="x0", cond=list(fac="1"), 
            main="... + s(fac) + s(fac, x0)", col='orange', ylim=c(7,22), rug=FALSE)
plot_smooth(gamm_int_slope, view="x0", cond=list(fac="2"), add=TRUE, col='red', xpd=TRUE)
plot_smooth(gamm_int_slope, view="x0", cond=list(fac="3"), add=TRUE, col='purple', xpd=TRUE)
plot_smooth(gamm_int_slope, view="x0", cond=list(fac="4"), add=TRUE, col='turquoise', xpd=TRUE)

Notez que les pentes aléatoires sont statique :

graphique des pentes aléatoires
plot(gamm_int_slope, select=3) 
# select=3 parce que la pente aléatoire se trouve sur la 3e ligne du tableau sommaire.

Enfin, nous allons examiner un modèle avec une surface lisse aléatoire.

Surface lisse aléatoire 
gamm_smooth <- gam(y ~ s(x0) + s(x0, fac, bs="fs", m=1), data=gam_data2)
summary(gamm_smooth)

Ici, si les pentes aléatoires variaient selon x0, nous auront vue des courbe variable pour chaque niveau :

graphiques des fonctions aléatoires
plot(gamm_smooth, select=1) 
# select=1 parce que le terme se trouve sur la 1e ligne du tableau sommaire.

Finalement, tous ces modèles mixes peuvent être compare en utilisant la fonction anova() pour trouver le meilleur modèle.

6. Autres distributions

Pour vous donner 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), l'exemple qui suit utilise un ensemble de données où une répartition binomiale sera nécessaire, y compris une modélisation d'une relation non linéaire. 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.

Loading the data
gam_data3 <- read.csv("other_dist.csv")
summary(gam_data3)
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  : Factor w/ 4 levels "f1","f2","f3",..: 1 1 1 1 1 1 1 1 1 1 ...

prop est la variable réponse, égal à la proportion de succès / (succès + échecs). Notez qu'il existe de nombreux cas où la proportion est égal à 1 ou 0 qui indique que les résultats ont toujours été des succès ou des échecs, respectivement, à ce moment mesuré durant l'expérience.
x1 est le temps écoulé depuis le début de l'expérience (variable explicative).
total représente le nombre de succès + échecs observé au moment x1i de l'expérience.
fac est un facteur qui code pour l'essai 1 à 4 de l'expérience (nous n'utiliserons pas cette variable dans cette section).

Commençons par la visualisation des données. Nous sommes intéressés par le nombre de succès par rapport aux échecs à mesure que x1 augmente. Étant donné qu'il y a des mesures répétées pour la valeur de x1 (essais 1 à 4, avec nombreuses observations par essai), nous pouvons d'abord présenter la proportion de succès en moyenne par boîte de temps (x1):

Visualisation des données
emptyPlot(range(gam_data3$x1), c(0,1), h=.5,
          main="Probability of successes", ylab="Probability",xlab="x1")
 
avg <- aggregate(prop ~ x1, data=gam_data3, mean, na.rm=TRUE)
lines(avg$x1, avg$prop, col="orange",lwd=2)

Notez comment la probabilité de succès augmente avec x1. D'après vous, est-ce que cette tendance est linéaire ou non linéaire? Nous allons tester cela en utilisant un GAM logistique (nous utilisons une distribution binomiale puisque la variable réponse représente des proportions).

GAM logistique
prop_model <- gam(prop~ s(x1), data=gam_data3, weights=total, family="binomial")
prop_summary <- summary(prop_model)
print(prop_summary$p.table)
print(prop_summary$s.table)
 
plot(prop_model)
             Estimate  Std. Error  z value   Pr(>|z|)
(Intercept)  1.173978  0.02709613  43.32641  0
        edf       Ref.df    Chi.sq     p-value
s(x1)   4.591542  5.615235  798.9407   2.027751e-164

Qu'est-ce que l'ordonnée représente dans ce modèle?

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

logit = log({N_{success}/N_{failures}})

- Si succès = échecs, le rapport est de 1 et le logit est 0 (log (1) = 0).
- Si les succès ont un nombre plus grand que les échecs, le ratio est supérieur à 1 et le logit a une valeur positive (par exemple, log(2) = 0,69).
- Si les succès ont un nombre plus petit que les échecs, le ratio est inférieur à 1 et le logit a une valeur négative (par exemple, log(0,5) = -0.69).

Donc, l'ordonnée est le logit, et indique s'il y a en moyenne plus de succès que d'échecs. Ici, l'estimé est positif ce qui signifie, qu'en moyenne, il n'y a plus de succès que d'échecs.

Qu'est-ce que le terme de lissage indique?

  • Ceci représente la façon dont les chances de succès vs échecs changent sur l'échelle de x1 (l'échelle du temps dans cet exemple). Donc, puisque l'edf > 1, la proportion de succès augmente plus rapidement au fil du temps (si par exemple, la réponse représente le nombre d'individus de l'espèce A vs l'espèce B et que nous augmentons la concentration des nutriments au fil du temps, ces résultats indiqueront que l'espèce A est de plus en plus observée alors que les concentrations de nutriments approchent de l'optimum de cette espèce au cours de l'expérience).

Visualiser la tendance au fil du temps

Enfin, nous allons voir les différentes façons de représenter ces relations graphiquement.

Tracer la tendance
par(mfrow=c(1,2))
plot(prop_model, select=1, scale=0, shade=TRUE)
abline(h=0)
 
plot_smooth(prop_model, view="x1",main="")
(diff <- find_difference(out$fv$fit, out$fv$CI, xVals=out$fv$x1))
addInterval(0, lowVals=diff$start, highVals = diff$end, col='red', lwd=2)
abline(v=c(diff$start, diff$end), lty=3, col='red')
text(mean(c(diff$start, diff$end)), 2.1, "sign. more \n success", col='red', font=3)

Quels renseignements ces graphiques nous apportent-ils vis à vis les succès et échecs ?

  • graphique de gauche: contribution (ou effet partiel si nous avions plus qu'une variable explicative) au fil du temps. La valeur logit augmente, donc les succès augmentent et les échecs diminuent.
  • graphique de droite: valeurs ajustées, ordonnée incluse (somme des effets si nous avions plus d'une variable explicative dans le modèle). Nous voyons ici que la valeur logit est estimée près de zéro au début de l'expérience ; cela signifie qu'il y a des quantités égales de succès et d'échecs. Peu à peu, les succès augmentent et à environ x1=400 il y a beaucoup plus de succès que d'échecs (l'effet est significativement différent de zéro). Nous avons également montré comment nous pouvons utiliser le graphique pour déterminer à quelle valeur de x1 cela se produit.

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

Graphique transformé
par(mfrow=c(1,1))
plot_smooth(prop_model, view="x1", main="",
            transform=plogis, ylim=c(0,1))
abline(h=.5, v=diff$start, col='red', lty=2)

Comme nous l'avons déjà vu avec le graphique précédent des valeurs logits, nous voyons qu'à approximativement x1=400 la proportion de succès augmente de façon significative au-dessus de 0,5.

7. Les GAMs en coulisse

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 d'une covariable, xi :

{y_i} = f(x_{i}) + {ε_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 choisissant une base, bi(x), définissant l'espace des fonctions dont f est un élément:

f(x) = sum{i=1}{q}{b_{i}(x)β_{i}}

Un exemple simple : une base polynomiale

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

b_{1}(x)=1, b_{2}(x)=x, b_{3}(x)=x^{2}, b_{4}(x)=x^{3}, b_{5}(x)=x^{4}

et f(x) devient:

f(x) = β_{1} + x_{i}β_{2} + x^{2}_{i}β_{3} + x^{3}_{i}β_{4}(x) + x^{4}_{i}β_{5}

et le modèle complet devient:

y_{i} = β_{1} + x_{i}β_{2} + x^{2}_{i}β_{3} + x^{3}_{i}β_{4}(x) + x^{4}_{i}β_{5} + ε_{i}

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

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

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

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 noeuds, 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 celle-ci fixée à 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}

il peut être modélisé en minimisant :

||y - XB||^{2} + {lambda}int{0}{1}{[f^{''}(x)]^{2}dx}

Quant lambda tend vers infty, le modèle devient linéaire. Pour la sélection du meilleur paramètre de lissage, lambda, on utilise une approche de validation croisée. 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 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 sum{i=1}{n}{(hat{f_{i}} - f_{i})^{2}}

Étant donné que f est inconnu, M est estimé en utilisant une technique de validation croisée généralisée qui laisse de côté, à chaque tour, une donnée et estime la capacité moyenne des modèles, construits sur les données restantes, de prédire la donnée qui a été mise de côté (pour plus de détails, consultez Wood 2006).

Illustration du principe de validation croisée

Dans le premier panneau, la courbe correspond à un ajustement faible par rapport aux données et ne fait pas mieux avec le point manquant. Dans le troisième panneau, la courbe ajuste le bruit aussi bien que le signal et la variabilité supplémentaire induite l'amène à prédire la donnée manquante plutôt mal. Dans le deuxième panneau, cependant, nous voyons que l'ajustement de la courbe du signal sous-jacent est très bien, le lissage passe à travers le bruit et la donnée manquante est raisonnablement bien prédite.

Note supplémentaire sur les degrés de liberté effectifs (edf)

Combien de degrés de liberté a un GAM ?

Au lieu de fournir la sortie de la validation croisée en termes de lambda (un paramètre qui détermine la complexité de l'ajustement), la fonction gam() utilise un terme appelé les degrés de liberté effectifs (edf), de manière cohérente à quantifier la complexité du modèle (pour justifier notre intuition que les degrés de liberté offrent une manière cohérente de quantifier la complexité du modèle). Parce que le nombre de paramètres libres des splines de lissage (tel que les GAMs) 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és. Donc, si le nombre de noeuds est arbitrairement réglé à k = 10, k-1 définit la limite supérieure des degrés de libertés associés à un terme de lissage. Ce nombre diminue alors que la pénalité lambda augmente jusqu'à ce que le meilleur modèle soit trouvé par validation croisée.

Références

Il existe beaucoup d'information sur les GAMs…

Simon Wood, l'auteur de la librairie mgcv, a un site très utile avec des conférences et des notes introductives sur la façon d'utiliser les GAMs :

Il a aussi écrit un livre, Generalized Additive Models: An Introduction with R, que nous avons utilisé comme référence pour cet atelier.

Le matériel de cet atelier a également été obtenu à partir des blogs et des tutoriels suivants :

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