Frage Ist es möglich, das Produkt von zwei glatten Begriffen in ein mgcv-Gam-Modell aufzunehmen?


Ich habe großen Erfolg mit Gam, um Saisonalität für Zeitreihen Daten zu modellieren. Mein neuestes Modell zeigt deutlich ein wöchentliches Muster zusätzlich zu saisonalen Änderungen. Während das wöchentliche Muster selbst über das Jahr sehr stabil ist, variiert seine Amplitude auch mit den Jahreszeiten. Im Idealfall möchte ich meine Daten wie folgt modellieren:

y ~ f(day in year) + g(day in year) * h(day in week)

woher f, g und h sind zyklische glatte Funktionen in mgcv

gam(
  y ~ s(day_in_year, k=52, bs='cc') 
  + s(day_in_year, k=52, bs='cc'):s(day_in_week, k=5, bs='cc')
  , knots=list(
    day_in_year=c(0, 356)
    , day_in_week=c(0,7)
  )
  , data = data
)

Leider funktioniert das nicht und wirft den Fehler zurück NA/NaN argument. Ich habe es versucht te(day_in_year, day_in_week, k=c(52, 5), bs='cc') das funktioniert, aber führt zu viele Freiheitsgrade ein, da das Modell Feiertage überlagert, die an bestimmten Wochentagen innerhalb der kurzen Anzahl von verfügbaren Jahren fallen.

Ist es möglich, ein Modell so zu spezifizieren, wie ich es versuche?


5
2018-05-12 22:45


Ursprung


Antworten:


Wow, eine ziemlich alte Frage!

In Bezug auf die Interaktion

Während das wöchentliche Muster selbst über das Jahr sehr stabil ist, variiert seine Amplitude auch mit den Jahreszeiten.

Die Verwendung von Tensorprodukt-Spline-Basis te ist der richtige Weg für die Interaktion, obwohl ein anständiger Konstruktor ist ti. Du hast das gesagt te gibt Ihnen viele Parameter zurück. Na sicher. Du hast k = 52 für den ersten Rand, und k = 5 für die Sekunde, dann enden Sie mit 52 * 5 - 1 Koeffizienten für diesen Tensorterm. Aber so entsteht eine Interaktion.

Beachten Sie, dass in mgcv GAM-Formel, : oder * gilt nur für die Interaktion zwischen parametrischen Termen. Die Interaktion zwischen Glättungen erfolgt über te oder ti.

Wenn das nicht das ist, worauf Sie hoffen, was erwarten Sie dann vom "Produkt"? Ein Hadamard-Produkt aus zwei marginalen Design-Matrizen? Was ist das für ein Gefühl? Ein Hadamard-Produkt benötigt übrigens zwei Matrizen gleicher Dimension. Ihre zwei Ränder haben jedoch nicht die gleiche Anzahl von Spalten.

Wenn Sie nicht verstehen, warum ich weiter über Matrizen spreche, dann müssen Sie Simons Buch im Jahr 2006 lesen. Obwohl die GAM-Schätzung erklärt, dass es jetzt ziemlich veraltet ist, wird die Konstruktion / Einrichtung von GAM (wie Designmatrizen und Strafmatrizen) in Kapitel erklärt 4 ändert sich auch nach einem Jahrzehnt überhaupt nicht.

OK, lassen Sie mich noch einen Hinweis geben. Das Reihe weise Kronecker Produkt verwendet für den Bau von te / ti Design-Matrix ist keine neue Erfindung.

Eine glatte Bezeichnung s(x) ist wie eine Faktorvariable g, als ob sie eine einzige Variable zu sein scheinen, sind sie konstruiert, um eine Entwurfsmatrix mit vielen Spalten zu sein. Zum g es ist eine Dummy-Matrix, während für f(x) es ist eine Basismatrix. So wird die Interaktion zwischen zwei glatten Funktionen auf die gleiche Weise konstruiert wie die Interaktion zwischen zwei Faktoren.

Wenn du einen Faktor hast g1 von 5 Ebenen und einem anderen Faktor g2 von 10 Ebenen, ihre marginale Design-Matrix (nach Kontrasten) haben 4 Spalten und 9 Spalten, dann die Interaktion g1:g2 würde 36 Spalten haben. Eine solche Design-Matrix ist nur das reihenweise Kronecker-Produkt der Design-Matrix von g1 und g2.


In Bezug auf Überanpassung

Wie Sie sagten, haben Sie nur wenige Jahre Daten, vielleicht 2 oder 3? In diesem Fall wurde Ihr Modell durch Verwendung von zu stark parametrisiert k = 52 zum day_in_year. Versuchen Sie es zum Beispiel zu reduzieren k = 30.

Wenn Überanpassung immer noch offensichtlich ist, hier sind ein paar Möglichkeiten, um es anzugehen.

Weg 1: Sie verwenden GCV zur Glättungsauswahl. Versuchen method = "REML". GCV neigt immer dazu, Daten zu überfrachten.

Weg 2: Bleiben Sie bei GCV und glätten Sie manuell die Glättungsparameter für eine schwerere Strafe. gamma Parameter von gam ist hier nützlich. Versuchen gamma = 1.5 beispielsweise.


Ein Tippfehler in deinem Code?

Die Knotenposition, sollte es sein day_in_year = c(0, 365)?


4
2018-05-17 07:41



Ihr Modell macht keinen großen Sinn, wenn Sie angeben, dass es

  1. ein saisonaler Effekt,
  2. ein Effekt der Wochentage,
  3. eine Interaktion, so dass der Wochentagseffekt mit dem Tag des Jahres variiert

Dies kann vollständig als ein Tensorprodukt glatt dargestellt werden. Das Modell, das Sie im Kommentar zur anderen Antwort erwähnen, finden Sie hier

y ~ f (Tag im Jahr) + g (Tag im Jahr) * h (Tag in der Woche)

ist nur eine Zerlegung des vollständigen Tensorprodukts ob was meinen Sie * als Haupteffekt + Interaktion. In diesem Fall ist das Modell, wie Sie es haben, nicht identifizierbar - Sie haben zweimal eine Funktion des Jahrestages. Und wenn du das Äquivalent von bedeutetest :, dann hat Ihr Modell nicht den Haupteffekt des Wochentages, was unerwünscht erscheint.

Ich passe Modelle dieser Form die ganze Zeit an (nur für Jahr und Tag des Jahres). Ich würde dies über folgendes erreichen:

gam(y ~ te(day_of_year, day_of_week, k = c(20, 6), bs = c("cc", "cc")),
    data = foo, method = "REML", knots = knots)

Sie könnten auch die Knoten-Definitionen optimieren. Ich neige dazu, Folgendes zu verwenden:

knots <- list(day_of_year = c(0.5, 366.5),
              day_of_week = c(0.5, 7.5)

Das wird keinen großen Unterschied machen, aber Sie setzen nur Randknoten nur ein wenig näher an die Daten.

Wenn Sie die Effekte zerlegen möchten, können Sie das Modell anpassen ti() glättet

gam(y ~ ti(day_of_year, bs = "cc", k = 12) + ti(day_of_week, bs = "cc", k = 6) + 
      te(day_of_year, day_of_week, k = c(12, 6), bs = c("cc", "cc")),
    data = foo, method = "REML", knots = knots)

Sie können die Werte von k um geeignete Werte in Verbindung mit Ihren Daten zu bestimmen gam.check().

Sie müssen dem Modell außerdem einen Begriff hinzufügen, um Ferien zu behandeln. Dies wäre ein parametrischer Term, der eine Anpassung anwendet, wenn der Tag ein Feiertag ist - also einen Faktor erstellen holiday und füge es dem Modell hinzu + holiday. Sie könnten an komplexere Modelle denken; Vielleicht ein Faktor Indexierung, wenn eine Woche einen Urlaub gekoppelt mit einem Faktor für die glatte hat day_of_week Komponente, so dass Sie ein wöchentliches Muster schätzen, wenn die Woche eine normale Woche ist, und ein zweites wöchentliches Muster, wenn die Woche einen Feiertag enthält.

Wenn Sie uns ein Beispiel / Diagramm der Daten zeigen, könnte ich weniger allgemein erweitern oder kommentieren.

Es ist nicht überraschend, dass die te() glatt passt Sie passt nicht gut in den Urlaub; Das Modell geht davon aus, dass der Effekt der Woche glatt und gleichmäßig mit einem reibungslosen Effekt des Tages des Jahres variiert. Feiertage sind keine weichen Abfahrten vom wöchentlichen Muster in der vorherigen oder nachfolgenden Woche. Das Urlaub Effekt wird nicht gut durch glatte Beziehungen modelliert und etwas anderes muss diesem Effekt Rechnung tragen.


3
2018-01-20 17:59