Différentes façons de modéliser les interactions entre les prédicteurs continus et catégoriels dans GAM

La question suivante s'appuie sur la discussion trouvée sur cette page . Étant donné une variable de réponse y, une variable explicative continue xet un facteur fac, il est possible de définir un modèle additif général (GAM) avec une interaction entre xet en facutilisant l'argument by=. Selon le fichier d'aide ?gam.models du package R mgcv, cela peut être accompli comme suit:

gam1 <- gam(y ~ fac +s(x, by = fac), ...)

@GavinSimpson suggère ici une approche différente:

gam2 <- gam(y ~ fac +s(x) +s(x, by = fac, m=1), ...)

Je joue avec un troisième modèle:

gam3 <- gam(y ~ s(x, by = fac), ...)

Mes principales questions sont: certains de ces modèles sont-ils tout simplement faux, ou sont-ils simplement différents? Dans ce dernier cas, quelles sont leurs différences? Sur la base de l'exemple que je vais discuter ci-dessous, je pense que je pourrais comprendre certaines de leurs différences, mais il me manque encore quelque chose.

À titre d'exemple, je vais utiliser un ensemble de données avec des spectres de couleurs pour les fleurs de deux espèces végétales différentes mesurées à différents endroits.

rm(list=ls())
# install.packages("RCurl")
library(RCurl) # allows accessing data from URL
df <- read.delim(text=getURL("https://raw.githubusercontent.com/marcoplebani85/datasets/master/flower_color_spectra.txt"))
library(mgcv)

Pour plus de clarté, chaque ligne de la figure ci-dessus représente le spectre de couleur moyen prévu pour chaque emplacement avec un GAM de forme distinct density~s(wl)basé sur des échantillons de ~ 10 fleurs. Les zones grises représentent 95% CI pour chaque GAM.

Mon objectif final est de modéliser l'effet (potentiellement interactif) Taxonet la longueur wld' onde sur la réflectance (appelés densitydans le code et l'ensemble de données) tout en tenant compte Localityd'un effet aléatoire dans un GAM à effets mixtes. Pour le moment, je n'ajouterai pas la partie effet mixte à mon assiette, qui est déjà suffisamment complète pour essayer de comprendre comment modéliser les interactions.

Je vais commencer par le plus simple des trois GAM interactifs:

gam.interaction0 <- gam(density ~ s(wl, by = Taxon), data = df) 
# common intercept, different slopes
plot(gam.interaction0, pages=1)

summary(gam.interaction0)

Produit:

Family: gaussian 
Link function: identity 

Formula:
density ~ s(wl, by = Taxon)

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  28.3490     0.1693   167.4   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
                      edf Ref.df     F p-value    
s(wl):TaxonSpeciesA 8.938  8.999 884.3  <2e-16 ***
s(wl):TaxonSpeciesB 8.838  8.992 325.5  <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.523   Deviance explained = 52.4%
GCV = 284.96  Scale est. = 284.42    n = 9918

La partie paramétrique est la même pour les deux espèces, mais des splines différentes sont adaptées pour chaque espèce. C'est un peu déroutant d'avoir une partie paramétrique dans le résumé des GAM, qui sont non paramétriques. @IsabellaGhement explique:

Si vous regardez les tracés des effets lisses estimés (lissages) correspondant à votre premier modèle, vous remarquerez qu'ils sont centrés sur zéro. Vous devez donc `` déplacer '' ces lissages vers le haut (si l'ordonnée à l'origine estimée est positive) ou vers le bas (si l'ordonnée à l'origine estimée est négative) pour obtenir les fonctions lisses que vous pensiez estimer. En d'autres termes, vous devez ajouter l'interception estimée aux lissages pour obtenir ce que vous voulez vraiment. Pour votre premier modèle, le «décalage» est supposé être le même pour les deux lissages.

Passons à autre chose:

gam.interaction1 <- gam(density ~ Taxon +s(wl, by = Taxon, m=1), data = df)
plot(gam.interaction1,pages=1)

summary(gam.interaction1)

Donne:

Family: gaussian 
Link function: identity 

Formula:
density ~ Taxon + s(wl, by = Taxon, m = 1)

Parametric coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)    40.3132     0.1482   272.0   <2e-16 ***
TaxonSpeciesB -26.0221     0.2186  -119.1   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
                      edf Ref.df    F p-value    
s(wl):TaxonSpeciesA 7.978      8 2390  <2e-16 ***
s(wl):TaxonSpeciesB 7.965      8  879  <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.803   Deviance explained = 80.3%
GCV = 117.89  Scale est. = 117.68    n = 9918

Désormais, chaque espèce a également sa propre estimation paramétrique.

Le modèle suivant est celui que j'ai du mal à comprendre:

gam.interaction2 <- gam(density ~ Taxon + s(wl) + s(wl, by = Taxon,  m=1), data = df)
plot(gam.interaction2, pages=1)

Je n'ai aucune idée claire de ce que ces graphiques représentent.

summary(gam.interaction2)

Donne:

Family: gaussian 
Link function: identity 

Formula:
density ~ Taxon + s(wl) + s(wl, by = Taxon, m = 1)

Parametric coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)    40.3132     0.1463   275.6   <2e-16 ***
TaxonSpeciesB -26.0221     0.2157  -120.6   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
                      edf Ref.df     F p-value    
s(wl)               8.940  8.994 30.06  <2e-16 ***
s(wl):TaxonSpeciesA 8.001  8.000 11.61  <2e-16 ***
s(wl):TaxonSpeciesB 8.001  8.000 19.59  <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.808   Deviance explained = 80.8%
GCV = 114.96  Scale est. = 114.65    n = 9918

La partie paramétrique de gam.interaction2est à peu près la même que pour gam.interaction1, mais maintenant il y a trois estimations pour les termes lisses, que je ne peux pas interpréter.

Merci d'avance à tous ceux qui prendront le temps de m'aider à comprendre les différences entre les trois modèles.

gam1et gam2vont bien; ce sont des modèles différents, bien qu'ils essaient de faire la même chose, c'est-à-dire des lissages spécifiques au groupe de modèles.

Le gam1formulaire

y ~ f + s(x, by = f)

pour ce faire, en estimant un lissage distinct pour chaque niveau de f(en supposant qu'il fs'agit d'un facteur standard), et en effet, un paramètre de lissage distinct est également estimé pour chaque lissage.

Le gam2formulaire

y ~ f + s(x) + s(x, by = f, m = 1)

atteint le même objectif que gam1(de modéliser la relation lisse entre xet ypour chaque niveau de f) mais il le fait en estimant un effet lisse global ou moyen de xsur y(le s(x)terme) plus un terme de différence lisse (le deuxième s(x, by = f, m = 1)terme). Comme la pénalité concerne ici la première dérivée ( m = 1) for this difference smoother, it is penalising departure from a flat line, which when added to the global or average smooth term (s (x) `) reflète un écart par rapport à l'effet global ou moyen.

gam3 forme

y ~ s(x, by = f)

est erronée quelle que soit sa capacité à s'adapter à une situation particulière. La raison pour laquelle je dis que c'est faux est que chaque lissage spécifié par la s(x, by = f)pièce est centré sur zéro en raison de la contrainte somme à zéro imposée pour l'identifiabilité du modèle. En tant que tel, il n'y a rien dans le modèle qui explique la moyenne de$Y$dans chacun des groupes définis par f. Il n'y a que la moyenne globale donnée par l'ordonnée à l'origine du modèle. Cela signifie que plus lisse, qui est centré sur zéro et dont la fonction de base plate a été supprimée de l'expansion de base de x(comme il est confondu avec l'ordonnée à l'origine du modèle) est maintenant responsable de la modélisation à la fois de la différence dans la moyenne de$Y$pour le groupe actuel et la moyenne globale (interception du modèle), plus l'effet lisse de xsur$Y$.

Cependant, aucun de ces modèles ne convient à vos données; en ignorant pour l'instant la mauvaise distribution de la réponse ( densityne peut pas être négatif et il y a un problème d'hétérogénéité qu'un non-gaussien familyrésoudrait ou résoudrait), vous n'avez pas pris en compte le regroupement par fleur ( SampleIDdans votre jeu de données).

Si votre objectif est de modéliser Taxondes courbes spécifiques, un modèle de formulaire serait un point de départ:

m1 <- gam(density ~ Taxon + s(wl, by = Taxon, k = 20) + s(SampleID, bs = 're'),
          data = df, method = 'REML')

où j'ai ajouté un effet aléatoire SampleIDet augmenté la taille de l'expansion de base pour les Taxonlissés spécifiques.

Ce modèle m1,, modélise les observations comme provenant soit d'un wleffet lisse selon l'espèce ( Taxon) d'où provient l'observation (le Taxonterme paramétrique définit simplement la moyenne densitypour chaque espèce et est nécessaire comme discuté ci-dessus), plus une interception aléatoire. Pris ensemble, les courbes pour les fleurs individuelles proviennent de versions décalées Taxondes courbes spécifiques, avec la quantité de décalage donnée par l'interception aléatoire. Ce modèle suppose que tous les individus ont la même forme de lissage que celle donnée par le lissé pour le particulier Taxondont la fleur individuelle provient.

Une autre version de ce modèle est la gam2forme ci-dessus mais avec un effet aléatoire supplémentaire

m2 <- gam(density ~ Taxon + s(wl) + s(wl, by = Taxon, m = 1) + s(SampleID, bs = 're'),
          data = df, method = 'REML')

Ce modèle correspond mieux mais je ne pense pas qu'il résout le problème du tout, voir ci-dessous. Je pense que cela suggère que la valeur par défaut kest potentiellement trop faible pour les Taxoncourbes spécifiques de ces modèles . Il y a encore beaucoup de variation douce résiduelle que nous ne modélisons pas si vous regardez les tracés de diagnostic.

Ce modèle est probablement trop restrictif pour vos données; certaines courbes de votre tracé des lissages individuels ne semblent pas être de simples versions décalées des Taxoncourbes moyennes. Un modèle plus complexe permettrait également des lissages spécifiques à chaque individu. Un tel modèle peut être estimé en utilisant la base d'interactionfs ou facteur-lisse . Nous voulons toujours Taxondes courbes spécifiques, mais nous voulons également avoir un lissage séparé pour chacune SampleID, mais contrairement aux bylissages, je suggère qu'au départ, vous voulez que toutes ces SampleIDcourbes spécifiques aient la même ondulation. Dans le même sens que l'interception aléatoire que nous avons incluse plus tôt, lafs base ajoute une interception aléatoire, mais inclut également une spline "aléatoire" (j'utilise les citations effrayantes comme dans une interprétation bayésienne du GAM, tous ces modèles ne sont que des variations sur des effets aléatoires).

Ce modèle est adapté à vos données comme

m3 <- gam(density ~ Taxon + s(wl, by = Taxon, k = 20) + s(wl, SampleID, bs = 'fs'), 
          data = df, method = 'REML')

Notez que j'ai augmenté kici, au cas où nous aurions besoin de plus de ondulation dans les lissages Taxonspécifiques. Nous avons encore besoin de l' Taxoneffet paramétrique pour les raisons expliquées ci-dessus.

Ce modèle prend beaucoup de temps pour s'adapter à un seul noyau avec gam()- bam()sera probablement mieux adapté à ce modèle car il y a un nombre relativement important d'effets aléatoires ici.

Si nous comparons ces modèles avec une version corrigée de la sélection des paramètres de lissage de l'AIC, nous voyons à quel point ce dernier modèle m3est nettement meilleur que les deux autres, même s'il utilise un ordre de grandeur plus de degrés de liberté.

> AIC(m1, m2, m3)
          df      AIC
m1  190.7045 67264.24
m2  192.2335 67099.28
m3 1672.7410 31474.80

Si nous regardons les lissages de ce modèle, nous avons une meilleure idée de la façon dont il ajuste les données:

(Notez que cela a été produit en draw(m3)utilisant la draw()fonction de mon package de gratia . Les couleurs dans le graphique en bas à gauche ne sont pas pertinentes et n'aident pas ici.)

SampleIDLa courbe ajustée de chacun est construite à partir de l'ordonnée à l'origine ou du terme paramétrique TaxonSpeciesBplus l'un des deux Taxonlissages spécifiques, en fonction de Taxonleur SampleIDappartenance, plus son propre SampleIDlissage spécifique.

Notez que tous ces modèles sont toujours erronés car ils ne tiennent pas compte de l'hétérogénéité; les modèles gamma ou Tweedie avec un lien de journal seraient mes choix pour aller plus loin. Quelque chose comme:

m4 <- gam(density ~ Taxon + s(wl, by = Taxon) + s(wl, SampleID, bs = 'fs'), 
          data = df, method = 'REML', family = tw())

Mais j'ai du mal avec ce modèle pour le moment, ce qui pourrait indiquer qu'il est trop complexe avec plusieurs lisses d' wlinclus.

Une autre forme consiste à utiliser l'approche factorielle ordonnée, qui effectue une décomposition de type ANOVA sur les lisses:

• Taxon le terme paramétrique est conservé
• s(wl)est un lissage qui représentera le niveau de référence
• s(wl, by = Taxon)aura une différence distincte lisse pour chaque autre niveau. Dans votre cas, vous n'en aurez qu'un.

Ce modèle est monté comme m3,

df <- transform(df, fTaxon = ordered(Taxon))
m3 <- gam(density ~ fTaxon + s(wl) + s(wl, by = fTaxon) +
            s(wl, SampleID, bs = 'fs'), 
          data = df, method = 'REML')

mais l'interprétation est différente; le premier s(wl)fera référence à TaxonAet le lissé impliqué par s(wl, by = fTaxon)sera une différence lisse entre le lissé pour TaxonAet celui de TaxonB.

C'est ce que Jacolien van Rij écrit dans sa page tutoriel:

La configuration de l'interaction dépend du type de prédicteur de regroupement:

• avec facteur incluent la différence d'interception: Group + s(Time, by=Group)
• avec facteur ordonné incluent la différence d'interception et la référence lisse: Group + s(Time) + s(Time, by=Group)
• avec prédicteur binaire incluent la référence lisse: s(Time) + s(Time, by=IsGroupChildren)

Les variables catégorielles doivent être spécifiées en tant que facteurs, facteurs ordonnés ou facteurs binaires avec les fonctions R appropriées. Pour comprendre comment interpréter les sorties et ce que chaque modèle peut et ne peut pas nous dire, consultez directement la page du didacticiel de Jacolien van Rij . Son didacticiel explique également comment adapter les GAM à effets mixtes. Pour comprendre le concept d'interactions dans le contexte des GAM, cette page de tutoriel de Peter Laurinec est également utile. Les deux pages fournissent de nombreuses informations supplémentaires pour exécuter correctement les GAM dans différents scénarios.