Verschiedene Methoden zur Modellierung von Interaktionen zwischen kontinuierlichen und kategorialen Prädiktoren in GAM

8

Die folgende Frage baut auf der Diskussion auf dieser Seite auf . Bei einer gegebenen Antwortvariablen y, einer kontinuierlichen erklärenden Variablen xund einem Faktor facist es möglich, ein allgemeines additives Modell (GAM) mit einer Interaktion zwischen xund unter facVerwendung des Arguments zu definieren by=. Gemäß der Hilfedatei ?gam.models im R-Paket mgcvkann dies wie folgt erreicht werden:

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

@ GavinSimpson schlägt hier einen anderen Ansatz vor:

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

Ich habe mit einem dritten Modell herumgespielt:

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

Meine Hauptfragen sind: Sind einige dieser Modelle einfach falsch oder sind sie einfach anders? Was sind im letzteren Fall ihre Unterschiede? Anhand des Beispiels, das ich unten diskutieren werde, denke ich, dass ich einige ihrer Unterschiede verstehen könnte, aber mir fehlt immer noch etwas.

Als Beispiel werde ich einen Datensatz mit Farbspektren für Blüten zweier verschiedener Pflanzenarten verwenden, die an verschiedenen Orten gemessen wurden.

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)

Aus Gründen der Klarheit repräsentiert jede Linie in der obigen Abbildung das mittlere Farbspektrum, das für jeden Ort mit einem separaten GAM der Form vorhergesagt wurde, density~s(wl)basierend auf Proben von ~ 10 Blumen. Die grauen Bereiche repräsentieren 95% CI für jedes GAM.

Mein letztes Ziel ist es, den (potenziell interaktiven) Effekt Taxonund die Wellenlänge wlauf das Reflexionsvermögen ( densityim Code und im Datensatz bezeichnet) zu modellieren und dabei Localityeinen zufälligen Effekt in einem GAM mit gemischten Effekten zu berücksichtigen . Im Moment werde ich den Mixed-Effect-Teil nicht zu meiner Platte hinzufügen, die bereits voll genug ist, um zu verstehen, wie Interaktionen modelliert werden.

Ich beginne mit dem einfachsten der drei interaktiven GAMs:

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

Geben Sie hier die Bildbeschreibung ein

summary(gam.interaction0)

Produziert:

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

Der parametrische Teil ist für beide Arten gleich, jedoch werden für jede Art unterschiedliche Splines angepasst. Es ist etwas verwirrend, einen parametrischen Teil in der Zusammenfassung der GAMs zu haben, die nicht parametrisch sind. @IsabellaGhement erklärt:

Wenn Sie sich die Diagramme der geschätzten Glättungseffekte (Glättungen) ansehen, die Ihrem ersten Modell entsprechen, werden Sie feststellen, dass sie um Null zentriert sind. Sie müssen diese Glättungen also nach oben (wenn der geschätzte Achsenabschnitt positiv ist) oder nach unten (wenn der geschätzte Achsenabschnitt negativ ist) verschieben, um die Glättungsfunktionen zu erhalten, von denen Sie dachten, dass Sie sie schätzen. Mit anderen Worten, Sie müssen den geschätzten Abschnitt zu den Glättungen hinzufügen, um das zu erreichen, was Sie wirklich wollen. Für Ihr erstes Modell wird angenommen, dass die Verschiebung für beide Glättungen gleich ist.

Weiter geht's:

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

Geben Sie hier die Bildbeschreibung ein

summary(gam.interaction1)

Gibt:

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

Jetzt hat jede Art auch ihre eigene parametrische Schätzung.

Das nächste Modell ist das, das ich nicht verstehen kann:

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

Geben Sie hier die Bildbeschreibung ein

Ich habe keine klare Vorstellung davon, was diese Grafiken darstellen.

summary(gam.interaction2)

Gibt:

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

Der parametrische Teil von gam.interaction2ist ungefähr der gleiche wie für gam.interaction1, aber jetzt gibt es drei Schätzungen für glatte Begriffe, die ich nicht interpretieren kann.

Vielen Dank im Voraus an alle, die sich die Zeit nehmen, mir zu helfen, die Unterschiede zwischen den drei Modellen zu verstehen.

Marco Plebani
quelle
2
Was für ein schöner Beitrag, Marco! Wenn Sie sich die Diagramme der geschätzten Glättungseffekte (Glättungen) ansehen, die Ihrem ersten Modell entsprechen, werden Sie feststellen, dass sie um Null zentriert sind. Sie müssen diese Glättungen also nach oben (wenn der geschätzte Achsenabschnitt positiv ist) oder nach unten (wenn der geschätzte Achsenabschnitt negativ ist) verschieben, um die Glättungsfunktionen zu erhalten, von denen Sie dachten, dass Sie sie schätzen. Mit anderen Worten, Sie müssen den geschätzten Abschnitt zu den Glättungen hinzufügen, um das zu erreichen, was Sie wirklich wollen. Für Ihr erstes Modell wird angenommen, dass die Verschiebung für beide Glättungen gleich ist.
Isabella Ghement
1
Wenn Sie Ihr Modell angeben, scheint es mir, dass Sie einen Haupteffekt für Taxon, einen Haupteffekt (glatt) für wl und eine (glatte) Interaktion zwischen Taxon und wl haben sollten. Der Link zu Gavin Simpsons Beitrag legt nahe, dass er auf diese Weise Modelle dieser Art erstellt. Er scheint auch den gleichen Wert von k für die glatten Effekte im Modell zu verwenden. Wenn Sie einen Interaktionsterm zwischen zwei Prädiktorvariablen einfügen, sollten Sie normalerweise auch die Haupteffekte für diese Variablen angeben.
Isabella Ghement
Daher würde ich Ihr erstes Modell verwerfen, da der Haupteffekt von Taxon weggelassen wird. Verwenden Sie einfach Gavins Vorschlag, um die wichtigsten Effekte und Interaktionseffekte zu erhalten, die Sie benötigen (und denken Sie daran, dass die vom Modell erzeugten Glättungen standardmäßig um 0 zentriert sind und je nach Intercept-Term (s) nach oben oder unten verschoben werden müssen.
Isabella Ghement
Hallo @IsabellaGhement und danke für dein Feedback. Wie würden Sie die Tatsache interpretieren, dass die Zusammenfassung (gam.interaction2) eine Signifikanzschätzung für s (wl) relativ zu jeder Art, aber auch eine für s (wl) ergibt, die mit keiner der beiden Arten verknüpft ist? Ist das die Auswirkung von wl auf die Glättungsfunktion von y (Dichte in meinem Fall) unabhängig von Taxon? Wird es einfach durch Anpassen der Dichte ~ s (wl) berechnet? Ich führe ein solches Modell aus und es schätzt einen parametrischen Koeffizienten. sehr nahe am Mittelwert des parametrischen Koeffizienten. der beiden Arten und des zugehörigen edf liegen sehr nahe an denen von s (wl), die in der Zusammenfassung (gam.interaction2) angegeben sind.
Marco Plebani
1
Kollegen und ich haben ein Papier in der Presse (Vorabdruck hier) , das viele Details zu diesen Themen enthält. Dies ist möglicherweise hilfreich, um sowohl die Auswahl an Modellen, die angepasst werden können, als auch die Auswahl zwischen diesen zu überprüfen. Für mich denke ich, dass alles, was Sie hier brauchen, gam1 plus etwas für den SampleIDEffekt ist und Sie etwas gegen das Problem der nicht konstanten Varianz tun müssen; Diese Daten scheinen aufgrund der Untergrenze nicht bedingt nach Gauß verteilt zu sein.
Gavin Simpson

Antworten:

1

gam1und gam2sind in Ordnung; Es handelt sich um verschiedene Modelle, obwohl sie versuchen, dasselbe zu tun, nämlich modellgruppenspezifische Glättungen.

Das gam1Formular

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

Dies geschieht durch Schätzen eines separaten Glätters für jede Stufe von f(unter der Annahme, dass dies fein Standardfaktor ist), und tatsächlich wird auch für jede Glättung ein separater Glättungsparameter geschätzt.

Das gam2Formular

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

erreicht das gleiche Ziel wie gam1(Modellierung der glatten Beziehung zwischen xund yfür jede Ebene von f), dies geschieht jedoch durch Schätzung eines globalen oder durchschnittlichen glatten Effekts von xon y(dem s(x)Term) plus eines glatten Differenzterms (dem zweiten s(x, by = f, m = 1)Term). Da die Strafe hier auf der ersten Ableitung liegt ( 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) `), spiegelt dies eine Abweichung vom globalen oder durchschnittlichen Effekt wider.

gam3 bilden

y ~ s(x, by = f)

s(x, by = f)YfxYxY

Keines dieser Modelle ist jedoch für Ihre Daten geeignet. Wenn Sie vorerst die falsche Verteilung für die Antwort ignorieren ( densitykann nicht negativ sein und es gibt ein Heterogenitätsproblem, das ein Nicht-Gaußscher familybeheben oder beheben würde), haben Sie die Gruppierung nach Blumen ( SampleIDin Ihrem Datensatz) nicht berücksichtigt .

Wenn Sie Taxonbestimmte Kurven modellieren möchten, ist ein Modell der Form ein Ausgangspunkt:

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

wo ich einen zufälligen Effekt für hinzugefügt SampleIDund die Größe der TaxonBasiserweiterung für die spezifischen Glättungen erhöht habe .

Dieses Modell m1modelliert die Beobachtungen so, dass sie entweder von einem glatten wlEffekt stammen, abhängig davon, von welcher Art ( Taxon) die Beobachtung stammt (der Taxonparametrische Term legt nur den Mittelwert densityfür jede Art fest und wird wie oben beschrieben benötigt), plus einem zufälligen Achsenabschnitt. Zusammengenommen ergeben sich die Kurven für einzelne Blumen aus verschobenen Versionen der Taxonspezifischen Kurven, wobei das Ausmaß der Verschiebung durch den zufälligen Achsenabschnitt gegeben ist. Dieses Modell geht davon aus, dass alle Individuen die gleiche Form von Glätte haben, wie sie durch die Glätte für die jeweilige TaxonBlume gegeben ist.

Eine andere Version dieses Modells ist die gam2Form von oben, jedoch mit einem zusätzlichen zufälligen Effekt

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

Dieses Modell passt besser, aber ich glaube nicht, dass es das Problem überhaupt löst, siehe unten. Eine Sache, die meiner Meinung nach nahe legt, ist, dass die Standardeinstellung kfür die Taxonspezifischen Kurven in diesen Modellen möglicherweise zu niedrig ist . Es gibt immer noch eine Menge glatter Restvariationen, die wir nicht modellieren, wenn Sie sich die Diagnosediagramme ansehen.

Dieses Modell ist höchstwahrscheinlich zu restriktiv für Ihre Daten. Einige der Kurven in Ihrem Diagramm der einzelnen Glättungen scheinen keine einfach verschobenen Versionen der TaxonDurchschnittskurven zu sein. Ein komplexeres Modell würde auch individuelle Glättungen ermöglichen. Ein solches Modell kann unter Verwendung der fsoder der faktorglatten Interaktionsbasis geschätzt werden . Wir wollen immer noch Taxonspezifische Kurven, aber wir wollen auch eine separate Glättung für jede SampleID, aber im Gegensatz zu den Glättungen bywürde ich vorschlagen, dass Sie anfangs möchten, dass alle diese SampleIDspezifischen Kurven die gleiche Wackeligkeit haben. Im gleichen Sinne wie der zufällige Abschnitt, den wir zuvor aufgenommen haben, ist derfs base fügt einen zufälligen Achsenabschnitt hinzu, enthält aber auch einen "zufälligen" Spline (ich verwende die Angstzitate wie in einer Bayes'schen Interpretation des GAM, alle diese Modelle sind nur Variationen zufälliger Effekte).

Dieses Modell ist für Ihre Daten als geeignet

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

Beachten Sie, dass ich hier zugenommen khabe, falls wir mehr Wackeligkeit in den spezifischen Glättungen benötigen Taxon. Wir brauchen den Taxonparametrischen Effekt aus den oben erläuterten Gründen immer noch.

Es dauert lange, bis dieses Modell auf einen einzelnen Kern passt gam()- bam()es ist höchstwahrscheinlich besser, dieses Modell anzupassen, da es hier eine relativ große Anzahl zufälliger Effekte gibt.

Wenn wir diese Modelle mit einer durch die Auswahl der Glättungsparameter korrigierten Version von AIC vergleichen, sehen wir, wie dramatisch besser dieses letztere Modell m3im Vergleich zu den beiden anderen Modellen ist , obwohl es eine Größenordnung mehr Freiheitsgrade verwendet

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

Wenn wir uns die Glättungen dieses Modells ansehen, erhalten wir eine bessere Vorstellung davon, wie es zu den Daten passt:

Geben Sie hier die Bildbeschreibung ein

(Beachten Sie, dass dies draw(m3)mit der draw()Funktion aus meinem Gratia- Paket erstellt wurde. Die Farben im Diagramm unten links sind irrelevant und helfen hier nicht weiter.)

Jede SampleIDangepasste Kurve wird entweder aus dem Achsenabschnitt oder dem parametrischen Term TaxonSpeciesBplus einer der beiden spezifischen Glättungen aufgebaut Taxon, je nachdem, zu welcher Taxonjeder SampleIDgehört, plus ihrer eigenen SampleIDspezifischen Glättung.

Beachten Sie, dass alle diese Modelle immer noch falsch sind, da sie die Heterogenität nicht berücksichtigen. Gamma- oder Tweedie-Modelle mit einem Protokolllink wären meine Wahl, um dies weiterzuentwickeln. Etwas wie:

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

Aber ich habe momentan Probleme mit dieser Modellanpassung, was darauf hindeuten könnte, dass sie zu komplex ist und mehrere Glättungen enthält wl.

Eine alternative Form ist die Verwendung des geordneten Faktoransatzes, der eine ANOVA-ähnliche Zerlegung der Glättungen durchführt:

  • Taxon Der parametrische Term bleibt erhalten
  • s(wl)ist eine Glättung, die den Referenzpegel darstellt
  • s(wl, by = Taxon)wird einen separaten Unterschied für jede andere Ebene glatt haben. In Ihrem Fall haben Sie nur eine davon.

Dieses Modell ist wie folgt ausgestattet 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')

aber die Interpretation ist anders; Das erste s(wl)bezieht sich auf TaxonAdas Glätten und das implizierte Glätten s(wl, by = fTaxon)ist ein glatter Unterschied zwischen dem Glätten für TaxonAund dem von TaxonB.

Gavin Simpson
quelle
Vielen Dank! Meine nächste Frage wäre gewesen: "Aber warum unterscheiden sich Zusammenfassungen, ob ein Faktor geordnet ist oder nicht?" aber du hast mich geschlagen, danke auch dafür. In meinem Datensatz ist jedes SampleIDein Spektrogramm von einer einzelnen Blume, jede von einer anderen Pflanze, daher denke ich nicht, dass SampleIDes als zufällig angegeben werden sollte (aber korrigiere mich, wenn ich falsch liege). Ich habe in der Tat ein Modell ähnlich Ihrem m3mit Taxonals geordnetem Faktor verwendet, aber + s(Locality, bs="re") + s(Locality, wl, bs="re")als zufällig angegeben. Ich werde mich mit den Fragen befassen, die Sie zur Verteilung von Residuen und zur Heteroskedastizität aufwerfen. Prost!
Marco Plebani
Ich würde immer noch SampleIDals zufällig angeben, dass die Daten einer einzelnen Blume wahrscheinlich in Beziehung stehen und mehr, wenn die gesamte Funktion, die sich auf die Blume bezieht, also in gewissem Sinne die Funktionen (Glättungen) zufällig sind. Möglicherweise benötigen Sie auch einen einfachen Zufallseffekt für eine Pflanze, wenn in der Studie mehrere Blüten pro Pflanze und mehrere Pflanzen pro Taxon vorhanden sind (verwenden Sie die bs = 're'"glatte", die ich zuvor in der Antwort erwähnt habe.
Gavin Simpson,
Als ich versuchte , passen m3mit family = Gamma(link = 'log')oder family = tw()ich war immer echte Probleme mit mgcv in der Lage, nicht gute Ausgangswerte und andere Fehler zu finden verursachen mgcv aus Mist, den ich nicht auf den Grund noch bekommen haben. Aus den von Ihnen angegebenen Daten ist ein Gaußsches Modell sicherlich nicht richtig. Ich habe einen Gaußschen mit Log-Link passend bekommen und es hat geholfen, aber es erfasst auch nicht die gesamte Heterogenität.
Gavin Simpson
0

Das schreibt Jacolien van Rij auf ihrer Tutorial-Seite:

Wie die Interaktion eingerichtet wird, hängt von der Art des Gruppierungsprädiktors ab:

  • mit Faktor umfassen Intercept-Differenz: Group + s(Time, by=Group)
  • mit geordnetem Faktor umfassen Intercept-Differenz und Referenzglättung: Group + s(Time) + s(Time, by=Group)
  • mit binärem Prädiktor enthalten Referenz glatt: s(Time) + s(Time, by=IsGroupChildren)

Kategoriale Variablen müssen als Faktoren, geordnete Faktoren oder binäre Faktoren mit den entsprechenden R-Funktionen angegeben werden. Um zu verstehen, wie die Ausgaben zu interpretieren sind und was jedes Modell uns sagen kann und was nicht, lesen Sie direkt die Tutorial-Seite von Jacolien van Rij . In ihrem Tutorial wird auch erklärt, wie man GAMs mit gemischten Effekten anpasst. Um das Konzept der Interaktionen im Kontext von GAMs zu verstehen, ist diese Tutorial-Seite von Peter Laurinec ebenfalls hilfreich. Beide Seiten bieten viele weitere Informationen, um GAMs in verschiedenen Szenarien korrekt auszuführen.

Marco Plebani
quelle