Binomial-Modell
Bisher haben wir hauptsächlich mit Samplingbasierter Inferenz, insbesondere MCMC gearbeitet. Nun wollen wir eine andere Art ansehen, wie man die Posteriori erhalten kann: Die Laplace-Approximation.
Uns liegen folgende Daten vor.
- Todesfälle nach Herzoperationen bei Babies
- Beobachtungen von 12 Krankenhäusern
- Annahme: Todesrate ist Krankenhausspezifisch
- Individuelle Raten der Krankenhäuser kommen aus gemeinamer Verteilung
Daten
Hierarchisches Modell
Datenmodell
Sei \(y_i\) die Anzahl der Todesfälle bei \(n_i\) Operationen in Krankenhaus \(i=1,\ldots,p=12\);
\[ \begin{aligned} y_i &\sim B(n_i, \pi_i); i=1,\ldots,p\\[3mm] \pi_i &= \frac{\exp(\theta_i)}{1-\exp(\theta_i)} \end{aligned} \]
Priori-Modell
Für den Effekt des Krankenhauses nehmen wir a priori eine Normalverteilung an:
\[ \theta_i|\tau \sim N\left(\mu,\tau^{-1}\right) \]
- Die Präzision \(\tau^{-1}\) ist dabei ein Parameter von Interesse: Wie stark weichen die Krankenhäuser voneinander ab?
- Wir wollen \(\tau^{-1}\) daher nicht vorgeben, sondern mitschätzen.
- Dafür brauchen wir eine (Hyper-)Priori.
Hyper-Priori
Wie bisher immer bei Präsizionen benutzen wir
\[ \tau \sim Ga(a,b) \]
z.B. mit \(a=b=0.001\).
Außerdem brauchen wir noch eine Hyper-Priori für \(\mu\), z.B. \(\mu\sim N(\mu_0,\sigma^2)\) mit hoher Varianz \(\sigma^2\) (relativ flache Priori).
Approximation
Versuchen wir nun die Posteriori zu approximieren. Wie beginnen ganz allgemein:
Hierarchische Bayes-Modelle
Sei allgemein ein HBM mit drei Stufen gegeben:
- \(y\sim F(\theta)\), \(\theta=(\theta_1,\ldots,\theta_p)\)
- \(\theta_j \sim N(\mu_j,\tau_j^{-1}),\, j=1, \ldots, p\)
- \(\tau_j \propto p(\tau_j),\, j=1,\ldots,p\)
Dabei ist \(F\) eine beliebige Datenverteilung.
Marginale Posteriori
Die marginale Posteriori von \(\theta|\tau\) ist
\[ p(\theta|y) = \int p(\theta|\tau, y) p(\tau| y) \; d\tau\\ \]
Hier brauchen wir für die Berechnung des Integrals die bedingte Posteriori von \(\theta\) und die marginale Posteriori von \(\tau\).
- Im Normalverteilungsmodell mit semi-konjugierter Priori kennen wir die bedingte Posteriori von \(\theta\) genau.
- Andernfalls müssten wir sie approximieren.
Die marginale Posteriori von \(\tau\) erhalten wir durch
\[ p(\tau|y) = \int p(\tau|\theta) p(\theta|y) d\theta \]
- Wir kennen diesen Ansatz aus dem Normalverteilungsmodell mit zwei unbekannten Parametern.
- In diesem hochdimensionalen Modell ist die Approximation allerdings etwas aufwendiger.
Laplace-Approximation
Betrachten wir erst mal allgemein dieses Problem: Gesucht ist das Integral
\[ \int_{-\infty}^{\infty} \exp(-nh(x))dx \]
Dabei sei \(h(x)\) eine konvexe, zweimal differenzierbare Funktion ist, die ihr Minimum an der Stelle \(x=\tilde{x}\) hat (und damit \(\exp(-nh(x))\) ihr Maximum).
Einfaches Beispiel: \(h(x)=(x-\mu)^2\) ist konvex und zweimal stetig differenzierbar. Dann ist \(\exp(-nh(x))\) proportional zur Dichte der Normalverteilung mit Erwartungswert \(\mu\).
Der Maximum-Likelihood-Schätzer maximiert \(\exp(-nh(x))\) und minimiert damit \(h(x)\).
Da \(\tilde{x}\) \(h(x)\) minimiert, gilt:
- \(h'(\tilde{x})=0\)
- \(h''(\tilde{x})>0\)
Wir können eine Taylorentwicklung von \(h(x)\) um \(\tilde{x}\) machen,, wobei wir nach dem zweiten Term abbrechen:
\[ \begin{aligned} h(x)& \approx h(\tilde{x})+h'(\tilde{x})(x-\tilde{x})+\frac12 h''(\tilde{x})(x-\tilde{x})^2\\ &= h(\tilde{x})+\frac12 h''(\tilde{x})(x-\tilde{x})^2 \end{aligned} \]
Setzen wir diese Taylor-Entwicklung ind die Exponentialfunktion ein, so erhalten wir
\[ \exp(h(x))\underset{\sim}{\propto}\exp\left(\frac12 h''(\tilde{x})(x-\tilde{x})^2\right) \]
- Das entspricht der Dichte der Normalverteilung.
- Alternativ kann man auch weitere Terme der Taylorreihe aufnehmen, z.B. kommt man mit dem dritten Term auf die sogenannte “schiefe Normalverteilung”.
Definition
Die Laplace-Approximation ist definiert als
\[ \begin{aligned} \int_{-\infty}^{\infty}\exp(-nh(x))dx &\approx \exp(-nh(\tilde{x})) \int_{-\infty}^{\infty} \exp\left(-\frac12 nh''(\tilde{x})(x-\tilde{x})^2\right)dx \\[3mm] &=\exp(-nh(\tilde{x}))\sqrt{\frac{2\pi}{nh''(\tilde{x})}} \end{aligned} \]
Der relative Fehler der Laplace-Approximation ist von der Ordnung \(O(\frac1n)\).
Laplace-Approximation des Posteriori-Erwartungswert
Berechnung des Posteriori-Erwartungswerts von \(g(\theta)\):
\[ \begin{aligned} E(g(\theta|x))&=\int_\Theta g(\theta) p(\theta|x) d\theta \\ &=\int_\Theta g(\theta) \frac{f(x|\theta)p(\theta)}{\int_\Theta f(x|\tilde{\theta}) p(\tilde{\theta}) d\tilde{\theta} }d\theta\\ &=\frac{\int_\Theta g(\theta)f(x|\theta)p(\theta)d\theta}{\int_\Theta f(x|\tilde{\theta})p(\tilde{\theta})d\tilde{\theta}} \end{aligned} \]
erfordert die Berechnung des Quotienten zweier ähnlicher Integrale.
Mit
\[ \begin{aligned} -nh(\theta)&=\log f(x|\theta)+\log p(\theta)\\ -nq(\theta)&=\log g(\theta)+\log f(x|\theta)+\log p(\theta) \end{aligned} \]
ist
\[ E(g(\theta)|x)=\frac{\int \exp(-nq(\theta))d\theta)}{\int \exp(-nh(\theta))d\theta)}. \]
Seien \(\hat{\theta}\) und \(\tilde{\theta}\) die Minimumstellen von \(h(\theta)\) und \(q(\theta)\). Dann ist
\[ E(g(\theta)|x)\approx \sqrt{\frac{h''(\hat{\theta})}{q''({\tilde{\theta}})}} \exp \left(-n(q(\tilde{\theta})-h(\hat{\theta})) \right). \]
Die Laplace-Approximation für \(E(g(\theta)|x)\) hat einen relativen Fehler von der Ordnung \(O(1/n^2)\).
Beispiel:Laplace-Approximation für das Poisson-Gamma-Modell
Simulierte Daten mit \(\lambda=25\) und Jeffreys-Priori:
| \(n\) | Posterior-Erwartungswert | Approximation | absoluter Fehler |
|---|---|---|---|
| 1 | 27.500 | 25.703 | 1.79745 |
| 2 | 28.750 | 27.801 | 0.94933 |
| 3 | 30.167 | 29.522 | 0.64487 |
| 5 | 26.300 | 25.909 | 0.39092 |
| 10 | 25.450 | 25.252 | 0.19763 |
| 25 | 25.700 | 25.620 | 0.07962 |
| 50 | 25.950 | 25.910 | 0.03991 |
| 100 | 26.175 | 26.155 | 0.01998 |
Für größere \(n\) funktioniert die Laplace-Approximation besser, weil die Binomial-Verteilung besser durch die Normalverteilung approximiert werden kann.
INLA
Für hierarchische Modelle kann man die Integrated Nested Laplace Approximation (INLA) benutzen.
Zur Erinnerung:
Sei allgemein ein HBM mit drei Stufen gegeben:
- \(y\sim F(\theta)\), \(\theta=(\theta_1,\ldots,\theta_p)\)
- \(\theta_j \sim N(\mu_j,\tau_j^{-1}),\, j=1, \ldots, p\)
- \(\tau_j \propto p(\tau_j),\, j=1,\ldots,p\)
Dabei ist \(F\) eine beliebige Datenverteilung.
Marginale Posteriori von \(\tau\)
Wir beginnen mit der Schätzung der marginale Posteriori von \(\tau\) und verwenden folgende Formel, die sich aus der Definition der bedingten Dichte ergibt:
\[ p(\tau|y) = \frac{p(\theta,\tau|y)}{p(\theta|\tau,y)} \]
\(p(\theta,\tau|y)\) ist bekannt bis auf die Normierungskonstante.
Für ein festes \(\tau\) können wir \(p(\theta|\tau,y)\) bestimmen:
- Für normalverteilten Daten kennen wir diese vollständig bedingte Posterior sogar (Normalverteilung).
- Andernfalls können wir sie mit der Laplace-Approximation approximieren (Notation: \(\tilde{p}(\theta|\tau,y)\)).
- Aus \(\tilde{p}(\theta|\tau,y)\) können wir einen Punktschätzer (Posteriori-Erwartungswert wie oben) für \(\theta^{*}(\tau)\) bestimmen und in die Formel einsetzen (damit die Approximation im vorherigen Schritt gut funktioniert).
Damit ergibt sich insgesamt als Approximation von \(p(\tau|y)\):
\[ \tilde{p}(\tau|y) = \left.\frac{p(\theta,\tau|y)}{\tilde{p}(\theta|\tau,y)}\right|_{\theta=\theta^*(\tau)} \]
Ähnlich wie beim CDF-Sampler können wir nun auf einem diskretisierten Gitter von \(\tau\) die Werte von \(\tilde{p}(\tau|y)\) bis auf die Normalisierungskonstante berechnen.
Daraus kann man die Normalisierungskonstante numerisch approxiomieren.
Oft haben wir mehrdimensionales \(\tau\), siehe z.B. Autofahrerbeispiel in den vorherigen Abschnitten.
Auch dann ist Diskretisierung möglich.
Rechenaufwand wächst jedoch exponential mit der Dimension, da Anzahl der Gitterpunkte exponential steigt.
Das ist der “Flaschenhals” der Methode.
Einige numerische Tricks können helfen, den Rechenaufwand zu minimieren:
- Gitter rund um grob geschätzes Maximum von \(p(\tau)\) aufbauen.
- Geschickt gewählte Gitterpunkte (Rotation der Dimensionen, möglichst Gitterpunkte mit hohem \(p(\tau)\) auswählen).
- Interpolation an Gitterpunkten statt expliziert Berechnung
Marginale Posteriori von \(\theta\)
Für die marginale Posteriori von \(\theta|y\) gilt nun:
\[ p(\theta|y) = \int p(\theta|\tau, y) p(\tau| y) d\tau \]
- Im vorherigen Schritt haben wir \(p(\tau| y)\) approximiert an festen \(\tau\).
- Dabei haben wir auch \(p(\theta|\tau, y)\) schon approximiert.
Damit können wir insgesamt die marginale Posteriori wie folgt maximieren:
\[ \tilde{p}(\theta|y)=\sum_\tau \tilde{p}(\theta|\tau,y)\tilde{p}(\tau|y) \]
R-INLA
Die INLA-Software kann die Approximation der Posteriori mit dem beschriebenen Verfahren durchführen. Wir benutzen das R-Paket INLA, welches die eigentliche Software an R anbindet.
Die Spezifikation des Modells ist ähnlich wie bei BayesX:
- Die Funktion f übernimmt dabei die Rolle der Modellierung der Priori.
- Hier ist hospital wieder die Nummer des Krankenhauses.
- “iid” steht für unabhängig und identisch verteilte Effekte pro Krankenhaus
- param=c(0.001,0.001) spezifiziert die Hyperparameter der Gamma-Priori (andere Prioris sind möglich)
Die Funktion inla führt den Algorithmus aus. family=“binomial” spezifiziert die Datenverteilung, Ntrials den Parameter \(n\) in der Binomialverteilung.
Ergebnisse
Wir erhalten aus dem Algorithmus die approximierte Posteriori
- des Intercepts
- der Effekte der einzelnen Krankenhäuser
- und der Präzision zwischen den Krankenhäusern
## Time used:
## Pre = 0.707, Running = 3.27, Post = 0.0447, Total = 4.02
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## (Intercept) -2.554 0.153 -2.878 -2.548 -2.266 -2.548 0
##
## Random effects:
## Name Model
## hospital IID model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant 0.975quant mode
## Precision for hospital 9.94 11.60 1.67 6.93 35.19 4.22
##
## Marginal log-Likelihood: -46.13
## is computed
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')