Uns liegen folgende Daten vor: Die Fälle von Lungenkrebs bei Männern in Westdeutschland zwischen 1975 und 1996 sowie die Bevölkerungsentwicklung in diesem Zeitraum.

Die Graphik zeigt die jährlichen Fälle in Tausend (\(y_t\)) und die Anzahl der erwarteten Fälle (\(e_t\)), errechnet aus \(b_t\) Bevölkerung im Jahr \(t\) errechnet durch:

\[ e_t=\frac{b_t\cdot\sum_ty_t}{\sum_tb_t} \]

Wir haben außerdem den geschätzten Tabakkonsum von Männern in Westdeutschland seit 1955. Tabakkonsum ist bekanntlich einer der Hauptgründe für Lungenkrebs. Die Graphik zeigt die Entwicklung seit 1955.

Wir wollen den Einfluß der Entwicklung des Tabakkonsums auf die Lungenkrebsfälle modellieren. Der Tabakkonsum “wirkt” dabei nicht sofort, sondern mit zeitlichem Abstand.

Da wir Zähldaten haben, gehen wir von einem Poisson-Modell aus. Für jedes Jahr setzen wir an:

\[ y_t \sim Po(\lambda_t e_t) \]

mit \(y_t\) Anzahl der Fälle im Jahr \(t\), \(e_t\) Anzahl der erwarteten Fälle im Jahr \(t\) und \(\lambda_t\) die unbekannte Rate im Jahr \(t\).

Nun wollen wir den Einfluß der Kovariablen Tabakkonsum auf die Fälle modellieren. Dabei folgen wir folgender Idee

• Im linearen Regressionmodell mit Normalverteilungsannahme haben wir festgestellt, das wir ein lineares Modell auf den Erwartungswert von \(y_t\) konstruiert haben.
• Der Erwartungswert der Poisson-Verteilung ist \(\lambda_te_t\). \(e_t\) ist gegeben, also modellieren wir \(\lambda_t\).
• Ein lineare Modell ist aber problematisch ist, da die Rate \(\lambda\) positiv ist.
• Idee: Wir transformieren \(\lambda\) mit dem Logarithmus und machen ein lineare Modell auf den Logarithmus:

\[ \log(\lambda_t) = \alpha+\beta x_t \]

Dabei ist

• \(\alpha\) der Intercept (Grundrate über alle Jahre) und
• \(\beta\) der lineare Einfluß der Kovariablen auf die Log-Rate.
• \(x_t\) ist der zeitverzögerte Tabakkonsum (z.B. 20 Jahre zuvor).

Priori

Im Poisson-Modell ist die Gamma-Verteilung die konjugierte Priori-Verteilung für \(\lambda\).

Hier allerdings haben wir \(\lambda\) transformiert. Wir übernehmen stattdessen die Prioris aus der linearen Regression:

• \(\alpha \sim N(m_\alpha,v_\alpha^2)\)
• \(\beta \sim N(m_\beta,v_\beta^2)\)

Datendichte

Leiten wir uns die Posteriori her. Für die einzelne Daten-Diche gilt:

\[ f(y_i|\lambda_i) = \frac{\lambda_i^{y_i}}{y_i!}\exp(-\lambda_i),\]

für die komplette Datendichte

\[ f(y|\lambda) = \prod_{t=1}^T\frac{\lambda_t^{y_t}}{y_t!}\exp(-\lambda_t) =\frac{\prod_{t=1}^T\lambda_t^{y_t}}{\prod_{t=1}^ny_t!}\exp\left(-\sum_{t=1}^n\lambda_t\right) \]

mit \(\lambda_t=\exp(\alpha+\beta x_t)\):

\[ f(y|\alpha,\beta) \propto \left(\prod_{t=1}^T\exp(\alpha+\beta x_t)^{y_t}\right)\exp\left(-\sum_{t=1}^T\exp(\alpha+\beta x_t)\right) =\\ = \exp \left(\alpha \sum_{t=1}^T y_t+\beta\sum_{t=1}^T x_ty_t-\sum_{t=1}^T\exp(\alpha+\beta x_t)\right) \]

Posteriori

Zusammen mit den Prioris ergibt sich die Posteriori

\[ \begin{align*} p(\alpha,\beta|y) \propto & \,\exp \left(\alpha \sum y_t+\beta\sum x_ty_t-\sum\exp(\alpha+\beta x_t)\right)\cdot\\ &\cdot\exp\left(-\frac{1}{2v_\alpha}(\alpha-m_\alpha)^2\right) \cdot\exp\left(-\frac{1}{2v_\beta}(\beta-m_\beta)^2\right)\\ =& \,\exp\left(-\frac{\alpha^2}{2v_\alpha}+\alpha\left(\frac{m_\alpha}{v_\alpha}+\sum y_t\right)\right)\\ &\cdot\exp\left(-\frac{\beta^2}{2v_\beta}+\beta\left(\frac{m_\beta}{v_\beta}+\sum x_ty_t\right)\right) \\&\cdot\exp\left(-\sum\exp(\alpha+\beta x_t)\right) \end{align*} \]

Man sieht, das bei diesem Modell Datendichte und Priori nicht im (semi-)konjugierten Sinne “zusammenpassen” . Der Gibbs-Sampler ist also nicht anwendbar (zumindest nicht ohne weiteres).

Wir kommen daher zu einem weiteren MCMC-Verfahren: Dem Metropolis-Algorithmus. Um aus einer mehrdimensionalen Verteilung \(f(\theta)\) zu ziehen, arbeitet er wie folgt:

1. Setzte Startwerte \(\theta^{(0)}\) fest
2. Setze \(k=1\)
3. Ziehe einen Vorschlag \(\theta^{*}\) aus einer symmetrischen Vorschlagsverteilung mit Dichte \(q(\theta^{*}|\theta^{(k-1)})\)
4. Berechne die Akzeptanzwahrscheinlichkeit \(\alpha=\min\left(1,\frac{f(\theta^{*})}{f(\theta^{(k-1)})}\right)\)
5. Ziehe \(u\) aus einer Gleichverteilung \(U[0,1]\)
6. Ist \(u\leq a\), dann akzeptiere den Vorschlag: \(\theta^{(k)}=\theta^{*}\)
7. Ist \(u>a\), dann verwerfe den Vorschlag und behalte den alten Wert: \(\theta^{(k)}=\theta^{(k-1)}\)
8. Setzte \(k=k+1\)
9. Wiederhole 3 bis 8 so oft wie nötig

Bemerkungen

• Symmetrische Vorschlagsverteilung heißt in diesem Fall

\[ q(\theta^{*}|\theta^{(k-1)})=q(\theta^{(k-1)}|\theta^{*}) \]

• Die einfachste Vorschlagsverteilung ist die Normalverteilung mit Erwartungswert \(\theta^{(k-1)}\), der sogenannte Random Walk.

• Ist die Dichte des Vorschlags größer als die Dichte des bisherigen Werts, dann wir der Vorschlag immer angenommen. Nur falls die Dichte des Vorschlags kleiner ist als die Dichte des bisherigen Werts, kann der Vorschlag verworfen werden.

• Wird der Vorschlag verworfen, dann wird der alte Wert wiederverwendet. Dass heißt, beim Metropolis-Algorithmus können mehrere Ziehungen hintereinander den selben Wert haben – auch wenn wir aus einer stetigen Verteilung ziehen.

• Bei der Berechnung der Akzeptanzwahrscheinlichkeit

\[ \alpha=\frac{f(\theta^{*})}{f(\theta^{(k-1)})} \]

fallen Konstanten weg, die nicht von \(\theta\) abhängen.

• Es reicht also, die Posteriori bis auf eine Proportionalitätskonstante zu kennen.
• Das ist praktisch, weil wir die Posteriori wirklich oft nur bis auf die Proportionalitätskonstante kennen.

Wie beim Gibbs-Sample erhalten wir auch hier abhängige Ziehungen aus der Posteriori von \(\alpha\) und \(\beta\). Hier die trace plots der ersten 5000 Iterationen, also gezogener Parameter geplottet gegen Iteration:

Die ersten circa 2000 Iterationen sind dabei wieder burn-in. Der Algorithmus braucht hier einige Zeit, bis er wirklich Ziehungen aus der Posteriori generiert. Die Startwerte wurden hier (absichtlich) ungünstig gewählt, so dass der Algorithmus einige Zeit braucht, bis er Paramter aus der richtigen Region zieht.

Die Länge des konkreten burn-in ist dabei – wie der ganze Algorithmus – zufällig. Später werden wir Möglichkeiten sehen, den burn-in besser zu bestimmen.

Wir lassen die ersten 5000 Iterationen weg und sehen uns die nächsten 5000 an.

Nun ist deutlich zu erkennen:

• Die Ziehungen sind von Iteration zu Iteration voneinander abhängig

• Manchmal wird über mehrere Iteration derselbe Wert gezogen

• Beides ist kein Problem, wir können aus den Ziehungen trotzdem die Posteriori-Dichte bzw. -Schätzer gewinnen

• Aber: der Algorithmus ist natürlich weniger effizient, als wenn wir unabhängige Ziehungen hätten.

• Der Metropolis-Algorithmus hängt von der verwendeten Vorschlagsverteilung (Proposal) ab.
• Hier benutzen wir die zweidimensionale Normalverteilung um den bisherigen Wert, eine alternative wäre zum Beispiel eine Gleichverteilung um den bisherigen Wert.
• Die Proposal-Verteilung hängt in der Regel von einer (oder mehrerer) Einstellungen ab, hier die Kovarianzmatrix der Normalverteilung.

Zum Vergleich die Trace-Plots der Iterationen 5001 bis 10000 mit besserer Einstellung der Kovarianzmatrix:

Man sieht an den trace plots:

• die Ziehungen scheinen “weniger abhängig” zu sein
• man erhält in weniger Iterationen Ziehungen, die über die komplette Posterioridichte verteilt sind; man sagt: der Algorithmus mischt besser (besseres Mixing)

ACF plots

Die Abhängigkeit zwischen den Ziehungen kann man sich durch die Auto-Korrelations-Funktion (auto correlation function, acf) visualisieren:

Dafür berechnet man für \(i=1,2,\ldots\) jeweils die Korrelation von Ziehungen, direkt hintereinander (\(i=1\)), von den jeweils zweitnächsten Ziehungen (\(i=2\)), usw. auf (\(i\) bezeichnet man als Lag):

Hier sehen wir (für \(\beta\)):

• bei Lag \(i=12\) liefgt die Autokorrelation noch über 0.5
• bei Lag \(i=25\) liegt die Autokorrelation bei ca. 0.25 liegt. Genauer: \(\beta^{(k)}\) und \(\beta^{(k+25)}\) haben eine Korrelation von \(\approx 0.248\).
• ab etwa Lag \(i=65\) ist die Autokorrelation fast gleich verschwunden, die Ziehungen mit Lag 65 sind fast voneinander unabhängig (die gestrichelte Linie zeigt die Grenzen eines Test aus Korrelation gleich Null).

Thinning

Geringere Autokorrelationen sprechen für einen effizienteren Algorithmus. Man könnte nun auf die Idee kommen, z.B. nur jede 65. Iteration zu benutzen. Man nennt dies Thinning, also Ausdünnen.

Als Beispiel hier der trace plot und die acf nach dem Ausdünnen mit dem Faktor 65, also nur jede 65. Ziehung wird verwendet:

• Die Ziehungen sind jetzt praktisch voneinander unabhängig.
• Aber Vorsicht: Wir verlieren dadurch jede Menge Ziehungen. Im Endeffekt werden unsere Schätzungen sehr viel schlechter.
• Thinning sollte daher nur verwendet werden, wenn man sowieso nicht alle Ziehungen speichern kann, z.B. in hochdimensionalen Modellen.

Schauen wir uns nun das Ergebnis unserer Analyse an.

Die folgende Tabelle zeigt

• Mittelwert,
• Median und
• 95%-Quantil der Ziehungen,

die uns

• Posteriori-Erwwartungwert,
• Posteriori-Median und
• symmetrisches 95%-Kredibilitätsintervall

approximieren.

	Mittelwert	Median	95%-Quantil
\(\alpha\)	9.828	9.827	[9.816, 9.840]
\(\beta\)	0.07965	0.07972	[0.07345, 0.08548]

• Posteriori-Modus und Highest-Posterior-Density-Intervall lassen sich nicht so einfach aus den Ziehungen gewinnen.
• Der Median als robustere Statistik wird oft gegenüber dem Mittelwert bevorzugt, letzterer lässt sich aber einfacher berechnen.

Für die vorliegenden Daten können wir also schliessen, dass eine Erhöhung des Tabakkonsums um den Wert 1 zu einer um (im Mittel) 0.797-fach höheren Log-Rate führt, also die Rate um \(\exp(0.797)\) den Faktor 2.219 multipliziert.

Die 0 ist im 95%-Kredibilitätsintervall nicht enthalten (auch nicht im 99%- oder 99,99-Kredibilitätsintervall), daher sind wir uns sehr sicher, dass es diesen Einfluß wirklich gibt.