Beispiel: Zeitreihe mit Saison-Effekt

Als Beispiel sehen wir und folgende Daten an:

Anzahl von getöteten oder schwer verletzten Autofahrern in England von Januar 1969 bis Dezember 1984

Beobachtungen

Wir sehen folgende Trends in den Daten:

• Es gibt einen Saisontrend, im Winter gibt es mehr Fälle als im Sommer
• Anfang 1983 sinkt die Anzahl der Fälle drastisch (zu diesem Zeitpunkt wurde die Gurtpflicht in Großbritannien eingeführt)
• Zwischen 1974 und 1983 bleibt die Anzahl der Fälle relativ konstant, davor schwankt sie stärker

Diese Trends wollen wir nun in Modell einbauen.

Datenmodell

Wir nehmen \(\sqrt{y}\) als normalverteilt an (für Fallzahlen wäre auch Poissonverteilung denkbar)

\[ \sqrt{y_i} \sim N(\mu_i, \sigma^2); i=1,\ldots,T=192 \]

Nun teilen wir den Erwartungswert \(\mu_i\) in verschiedene additive Effekte auf:

\[ \mu_i = \alpha + \beta x_i + \gamma_i + \delta_i \]

mit \(i\) Monat, \(x_i\) Dummyvariable Gurtpflicht ja/nein.

Dabei ist

• \(\alpha\) der Intercept
• \(\beta\) der Einfluß der Gurtpflicht
• \(\gamma\) der allgemeine Zeittrend
• \(\delta\) der Saisontrend

Wir wollen die zeitliche Entwicklung modellieren.

• Eine Lineare Entwicklung wäre vermutlich falsch
• Eine Parametrische Modellierung wäre eventuell möglich, aber welches Modell würden wir wählen?

Wir nehmen an, dass (wenn wir saisonalen Trend und Effekt der Gurtpflicht rausrechnen) aufeinander folgende Monate ähnlich viele Fälle haben.

• Idee: Dummykodierung, ein Parameter pro Monat
• \(\gamma_i\) ist der Effekt des Monats \(i\)

Der Monat \(i\) hat ähnlich viele Fälle wie der Monat \(i-1\), das ist unsere Priori-Annahme. Bringen wir das in Form einer Verteilung, genauer einer Normalverteilung:

\[ \gamma_i|\gamma_{i-1},\tau_c \sim N(\gamma_{i-1},\tau_c^{-1}); i=2,\ldots,T \]

• Im Monat \(i\) erwarten wir also so viele Fälle wie im Monat \(i-1\), mit einer Varianz \(\tau_c^{-1}\).
• Für den ersten Monat nehmen wir eine nicht-informative Priori: \(p(\gamma_1)\propto const.\).
• Man nennt diesen Ansatz auch Random Walk (zufälliger Gang): Jeder Schritt weicht, unabhängig von den vorherigen Schritten, vom momentanen Punkt zufällig nach oben oder unten ab.
• Man spricht auch von Random Walk erster Ordnung (RW1), höhere Ordnung beziehen mehrere vergangene Zeitpunkte mit ein.

Nun haben wir einen Parameter in der Priori: \(\tau_c\), die inverse Varianz (Präzision) in unserer Priori-Annahme. Um so größer die Präzision (damit umso kleiner die Varianz), umso ähnlicher sind sich die Monate. Ist die Präzision klein, dann ist die Varianz in unserem zeitlichen Trend größer.

Wie legen wir diesen Prioriparameter \(\tau_c\) fest?

• Wir wollen den Prioriparameter nicht festlegen, wir wollen ihn als unbekannt behandeln.
• Unbekannte Parameter haben im Bayesianischen eine Verteilung.
• Wir legen also eine Priori-Verteilung fest und erhalten nach Beobachtung eine Posteriori-Verteilung für \(\tau_c\).
• \(\tau_c\) ist eine inverse Varianz im einer Normalverteilung. Dafür haben wir bisher immer die Gamma-Verteilung benutzt:

\[ \tau_c\sim Ga(a_c,b_c) \]

Hierarchisches Modell

Wir haben nun ein sogenanntes Hierarchisches Bayes-Modell (HBM). Ein hierarchisches Modell hat mehrere Schichten, in unserem Beispiel (und meistens) drei:

• Level 1: Datenmodell, Definition der Likelihood
• Level 2: Priori-Modell der unbekannten Parameter
• Level 3: (Hyper-)Prioris der Prioriparameter in Level 2

Hyperpriori

Auch hier wollen wir die Präzision \(\tau_d\), also die inverse Varianz des Saisoneffekts, nicht festlegen, sondern mitschätzen. Wir geben wieder eine Hyperpriori vor:

\[ \tau_d\sim Ga(a_d,b_d) \]

Priori-Modell - Intercept und Kovariableneffekt

• Es fehlt noch die Priori-Information für den Intercept \(\alpha\) und den Kovariableneffekt \(\beta\).
• Wie im Regressionsmodell haben wir keine Priori-Information

\[ \begin{aligned} p(\alpha)& \propto \text{const.}\\ p(\beta)& \propto \text{const.} \end{aligned} \]

Priori auf Varianz

Schliesslich brauchen wir noch eine Priori für die Varianz in der Normalverteilung von \(\sqrt{y_i}\): \(\sigma^2\sim IG(a_s,b_s)\)

Insgesamt erhalten wir:

\[ \begin{aligned} p(\alpha,\beta,\gamma,\delta,\tau_c,\tau_d,\sigma^2|y^*) \propto &\, f(y^*|\alpha,\beta,\gamma,\delta,\sigma^2)\\ &\cdot p(\alpha)\cdot p(\beta)\\ &\cdot p(\gamma|\tau_c)p(\tau_c)\\ &\cdot p(\delta|\tau_d)p(\tau_d)\\ &\cdot p(\sigma^2) \end{aligned} \]

Full conditional der Priori-Parameter

Für einen MCMC-Algorithmus brauchen wir wieder die Full Conditionals. Sehen wir uns die Full Conditional von \(\tau_c\) an:

\[ \begin{aligned} p(\tau_c|\cdot) \propto &\, p(\alpha,\beta,\gamma,\delta,\tau_c,\tau_d,\sigma^2|y)\\[2mm] \propto &\, f(y^*|\alpha,\beta,\gamma,\delta,\sigma^2)\\ &\cdot p(\alpha)\cdot p(\beta)\\ &\cdot p(\gamma|\tau_c)p(\tau_c)\\ &\cdot p(\delta|\tau_d)p(\tau_d)\\ &\cdot p(\sigma^2)\\[2mm] \propto &\, p(\gamma|\tau_c)p(\tau_c) \end{aligned} \]

• Die vollständig bedingte Posteriori von \(\tau_c\) hängt also gar nicht von den Daten ab.
• Die marginale Posteriori von \(\tau_c\) hängt – über \(\gamma\) – natürlich schon von den Daten ab.
• Jedes hierarchische Modell ist so aufgebaut, dass jede Ebene nur von der Ebene darüber und der Ebene darunter abhängt.
• Noch mehr: \(\tau_c\) hängt a posteriori nur vom zugehörigen \(\tau\) ab, nicht von den anderen Parametern.

Hier kann man zeigen, dass die Gamma-Priori von \(\tau_c\) die semi-konjugierte Priori ist, also auch die Full Conditional eine Gamma-Verteilung ist.

Die komplette Herleitung der Full Conditionals und die Implementierung des Modells finden Sie hier.

Kovariableneffekt

Folgende Graphik zeigt die Posteriori des Effekts der Gurtpflicht:

• Der Posteriori-Median des Effekts der Gurtpflicht liegt bei -3.078, der Posteriori-Erwartungswert liegt bei -3.076.
• Wir erwarten im Mittel einen Effekt der Gurtpflicht eine Reduzierung der Anzahl der Fälle von 9.720. Das ist der der Posteriori-Erwartungswert von \(\beta^2\).

Zeiteffekt

Pro Monat (\(T=192\)) haben wir einen Zeiteffekt geschätzt. Wie können hier nicht alle Posterioridichten anschauen, stattdessen zeigen wir Median und 95%-Kredibilitätsintervall pro Monat:

Wir sehen den Anstieg bis etwa 1973, den starken Abfall bis etwa 1976, danach nur mehr ein leichter Trend nach unten. Der starke Abfall der Daten 1983 ist im Zeittrend nicht zu finden, sondern wird durch beta erklärt.

Wir können uns auch die Differenzen der Zeiteffekte ansehen. Zur Erinnerung, unsere Priori war:

\[ \gamma_t\sim N(\gamma_{t-1},\tau_c^{-1}) \]

oder auch

\[ \gamma_t-\gamma_{t-1}\sim N(0,\tau_c^{-1}) \]

Die Differenzen haben also a priori den Erwartungswert 0. Sehen wir uns Median und 95%-Kredibilitätsintervall der Differenzen an:

• Generall schwanken die Differenzen um (den Erwartungswert) 0.
• Das heißt, der Zeiteffekt bleibt relativ konstant.
• Rasche Veränderungen der Daten (z.B. Herbst 1973 und Herbst 1974) zeigen sich in der Posteriori der Differenzen.

Saisoneffekt

Auch für den Saisoneffekt hatten wir einen Effekt pro Monat:

• Der Saisontrend ist deutlich erkennbar.
• November und Dezember stellen sich als gefährlichste Monate heraus (etwa +5 im Mittel), April als ungefährlichster (bis zu -3).
• Hier kann der Saisoneffekt von Jahr zu Jahr sich ändern.
• A priori bleibt die Summe von 12 Monaten im Erwartungswert gleich Null.

Prioriparameter

Die beiden Prioriparameter \(\tau_c\) und \(\tau_d\) haben folgende Posterioris:

Diese Prioriparameter lassen sich wie folgt definieren:

• \(\tau_c\) entspricht der inversen Varianz der \(\gamma\)-Differenzen

• \(\tau_d\) entspricht der inversen Varianz der Summe von zwölf aufeinanderfolgenden Monaten für \(\delta\)

• Zum Beispiel werden für höheres \(\tau_c\) die Differrenzen immer kleiner, damit die Zeitreihe der \(\gamma\)-Effekte immer mehr zu einer konstanten Gerade.

• Kleineres \(\tau_c\) spricht für höhere Varianz der Differenzen und damit stärkeres Schwanken des Zeiteffekts.

• Wichtig: wir haben nicht einen \(\tau_c\)-Wert, sondern eine Posteriori-Verteilung von \(\tau_c\). Die Posteriori-Verteilung von \(\gamma\) berücksichtigt diese Posteriori von \(\tau_c\) (und natürlich umgekehrt), d.h. die Unsicherheit über \(\tau_c\) geht auch in \(p(\gamma|y)\) ein.

Glatte Zeitreihe

• Zuletzt können wir noch alle Effekte aufaddieren und das Ergebnis quadrieren.
• Wir erhalten eine Schätzung der Fälle.
• Dies muss pro Ziehung passieren! Dann erhalten wir sogar eine Posteriori für jedes \(y_t\).

Wir haben bei diesem hierarchischen Modell folgende Paramter

• \(\alpha\) und \(\beta\),
• jeweils \(T=172\) \(\gamma_t\) und \(\delta_t\)-Parameter,
• \(\tau_c\) und \(\tau_d\) sowie
• \(\sigma^2\).

• Insgesamt also \(2+2*192+2+1=389\) Parameter!
• Zweifellos ein hochdimensionales Modell.
• Uns liegen \(T=172\) Beobachtungen vor! Sehr viel weniger als Parameter
• Die Schätzung kann hier nur funktionieren, weil wir zusätzliche Priori-Information ins Modell stecken.

• Die Modellierung der Daten wurde hier in die Priori verlagert, zum einen durch einen relativ “glatten” Zeiteffekt, zum anderen durch einen Saisoneffekt.
• Hierarchische Modellierung ist sehr flexibel in der Modellierung
• Dadurch, dass die einzelnen Ebenen bzw. sogar die einzelnen Parameter nur von den direkt “benachbarten” Parametern abhängen, wird die Inferenz nicht sehr viel komplexer.