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

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)

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)

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)

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.

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)$Y$fx$Y$x$Y$

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:

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

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.