Skip to Tutorial Content

Empirischer Bayes

In einem früheren Abschnitt hatten wir folgende Daten:

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

und hatten folgendes Modell angenommen:

\[ \begin{aligned} \sqrt{y_i} & \sim N(\mu_i, \sigma^2);\, i=1,\ldots,T=192 \\[2mm] \mu_i & = \alpha + \beta x_i + \gamma_i + \delta_i \end{aligned} \]

  • Für die Effekte \(\theta=(\alpha, \beta, \gamma_i, \delta_i)\) hatten wir jeweils Normalverteilungspriori angenommen.
  • Für den MCMC-Algorithmus ergab sich ein Gibbs-Sampler, da dies die semi-konjugierte Prioris sind, z.B.

\[ \theta|\tau \sim N(0,\tau) \Rightarrow \theta|\cdot\sim N(\tilde{\mu},\tilde{\tau}) \]

  • Das gilt nicht nur für die einzelnen Parameter, auch \(\theta\) insgesamt ist multivariat normalverteilt.
  • Würden wir \(\tau\) kennen, würden wir die Posteriori komplett kennen!

BayesX

Die Software BayesX schätzt für bestimmte Arten von Modellen realisiert mittels MCMC oder mittels empirischer Bayes-Ansätze.

BayesX is a software tool for estimating structured additive regression models. Structured additive regression embraces several well-known regression models such as generalized additive models (GAM), generalized additive mixed models (GAMM), generalized geoadditive mixed models (GGAMM), dynamic models, varying coefficient models, and geographically weighted regression within a unifying framework…

BayesX Logo

BayesX

Wir verwenden Daten und Modell aus dem vorherigen Abschnitt und das Paket R2BayesX.

library(R2BayesX)

Zuerst müssen wir das Modell spezifizieren:

formula.bx  = sqrt(y) ~ belt + sx(trend, bs="rw1", a=1, b=0.0005) +  sx(seasonal, bs="season", a=1, b=0.1)

In der Variable formula.bx speichern wir das Modell ab:

  • sqrt(y) soll modelliert werden

  • belt ist eine Kovariable

  • trend und seasonal sind einfach die laufenden Monatszahlen

  • sx spezifiziert Priori-Modelle:

    • bs=“rw1” bezeichnet unser RW1-Priori
    • bs=“season” entspricht unserer Saison-Priori
    • Die Hyperpriori-Parameter werden mit übergeben.

MCMC

Die Funktion bayesx führt den Algorithmus durch. Wir wählen hier method=“MCMC” und müssen noch die Verteilung der Daten (“gaussian”) und Hyper-Parameter für die Varianz der Daten sowie die Anzahl der Iterationen (und des burn-in) angeben.

mcmc.bx <- bayesx(formula = formula.bx, data=Drivers, 
      control=bayesx.control(family="gaussian", method="MCMC", hyp.prior = c(4,4), 
      iterations=10000L, burnin=5000L))

Die Ergebnisse des Algorithmus lassen sich mit plot bzw. summary ausgeben.

plot(mcmc.bx$effects$`sx(trend)`)

plot(mcmc.bx$effects$`sx(seasonal)`)

summary(mcmc.bx)
## Call:
## bayesx(formula = formula.bx, data = Drivers, control = bayesx.control(family = "gaussian", 
##     method = "MCMC", hyp.prior = c(4, 4), iterations = 10000L, 
##     burnin = 5000L))
##  
## Fixed effects estimation results:
## 
## Parametric coefficients:
##                Mean      Sd    2.5%     50%   97.5%
## (Intercept) 41.2508  0.1727 40.9177 41.2432 41.6039
## belt        -4.4804  1.0568 -6.6987 -4.3893 -2.5518
## 
## Smooth terms variances:
##                Mean     Sd   2.5%    50%  97.5%    Min    Max
## sx(seasonal) 0.2379 0.0865 0.1141 0.2232 0.4512 0.0925 0.8240
## sx(trend)    0.1976 0.0934 0.0771 0.1763 0.4212 0.0494 1.0161
##  
## Scale estimate:
##          Mean     Sd   2.5%    50%  97.5%
## Sigma2 2.8205 0.3351 2.2539 2.7922 3.5351
##  
## N = 192  burnin = 5000  method = MCMC  family = gaussian  
## iterations = 10000  step = 10

BayesX mit REML

Nun das selbe Modell mit dem empirischen Bayes- bzw. REML-Ansatz.

Die Definition des Modells ist identisch, außer dass keine Hyperprioris spezifiziert werden!

Mit method=“REML” erhalten wir den REML-Schätzer:

##    user  system elapsed 
##   3.568   0.044   3.760

## Call:
## bayesx(formula = formula.bx, data = Drivers, control = bayesx.control(family = "gaussian", 
##     method = "REML"))
##  
## Fixed effects estimation results:
## 
## Parametric coefficients:
##             Estimate Std. Error  t value  Pr(>|t|)    
## (Intercept)  41.2769     0.1665 247.8527 < 2.2e-16 ***
## belt         -4.6514     1.0544  -4.4114 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Smooth terms:
##              Variance Smooth Par.       df Stopped
## sx(seasonal)   0.0193    117.0120  14.2552       0
## sx(trend)      0.2541      8.8647  19.0758       0
##  
## Scale estimate: 2.2527 
##  
## N = 192  df = 35.331  AIC = 383.257  BIC = 498.347  
## GCV = 2.76067  logLik = -156.2975  method = REML  family = gaussian

  • Für die Präzisionsparameter (Smooth terms) erhalten wir nur Punktschätzer, keine Posteriori
  • Die Intervallschätzer der Effekte \((\theta)\) sind beim empirischen Bayes-Ansatz tendenziell kleiner (dazu gleich mehr)
  • Davon abgesehen unterschieden sich die Ergebnisse nur in Details.

Model Averaging

Warum sind die Intervallschätzer beim Empirischen Bayes-Ansatz kleiner?

  • Der empirische Bayes-Ansatz berücksichtigt die Unsicherheit über die Präzisionsparameter \(\tau\) nicht.
  • Die Unsicherheit wird durch die Posteriori von \(\tau\) ausgedrückt.
  • Nur im “vollen Bayes-Ansatz” (z.B. mit MCMC) erhalten wir die Posteriori von \(\theta\) und \(\tau\) und berücksichtigen in der Schätzung von \(\theta\) auch die Unsicherheit über \(\tau\).

  • Die Schätzung der \(\tau\) ist in der Regel nicht besonder interessant.
  • Uns interessiert die marginale Posterioriverteilung von \(\theta\). Für diese gilt:

\[ p(\theta|y)=\int p(\theta,\tau|y) d\tau = \int p(\theta|\tau,y)p(\tau|y) d\tau \]

  • Im Gegensatz zum Empirischen Bayes-Ansatz erhalten wir beim “vollen Bayes-Ansatz” also nicht das Ergebnis für einen \(\tau\)-Wert, sondern die Mischung von verschiedenen Modellen, gewichtet mit der marginalen Posteriori-Verteilung von \(\tau\).
  • Man spricht hier auch von Model Averaging.

Weiter

BayesX