Um Bayes-Modelle umzusetzen, gibt es verschiedene Software-Pakete. Dabei gibt es ziemlich universal einsetbare wie BUGS, JAGS oder STAN oder auch spezialisierte wie BayesX oder INLA.

Ansätze

All diese Pakete basieren auf unterschiedlichen Ansätzen von MCMC oder auch Approximation:

• BUGS (Bayesian inference Using Gibbs Sampling) benutzt nur Gibbs-Sampler, wobei die full conditionals numerisch approximiert werden
• STAN nutzt eine moderne MCMC-Variante namens Hamiltonian Monte Carlo (oder alternativ diverse Approxmiations-Verfahren)
• INLA (Integrated Nested Laplace Approximation) verbindet mehrere Approximationsverfahren (insbesondere Laplace-Approximation) und ist für additive Regressionsmodelle geeignet.

BUGS und JAGS

• BUGS startete 1989 als WinBUGS. Seit 2005 wurde das Projekt langsam in OpenBUGS überführt.

• JAGS (Just another Gibbs sampler) ist eine Re-Implementierung von BUGS (geschrieben in C++).

• Die Modellierungssprache in BUGS und JAGS ist praktisch identisch.

• In R gibt es mehrere Pakete für OpenBUGS als auch für JAGS, wir benutzen hier rjags:

library(rjags)

## Loading required package: coda

## Linked to JAGS 4.3.0

## Loaded modules: basemod,bugs

JAGS arbeitet mit einer einfach aufgebauten Modellierungssprache.

Erinnern wir uns an unser Poisson-Regressions-Modell. Unser Modell lautete:

\[ \begin{eqnarray} y_t &\sim& Po(\lambda_t e_t)\\ \log(\lambda_t) &=& \alpha+\beta x_t\\ \alpha &\sim& N(m_\alpha,v_\alpha^2)\\ \beta &\sim& N(m_\beta,v_\beta^2) \end{eqnarray} \]

Diese mathematische Modellbeschreibung setzen wir jetzt in die JAGS-Sprache um:

model <- "model{
    for (t in 1:T) {
        y[t] ~ dpois(lambda[t]*e[t])
        log(lambda[t]) <- alpha + beta*x[t]
    }
    # Prioris
    alpha ~ dnorm(0, 0.001)
    beta ~ dnorm(0, 0.001)
}"

Die Modellbeschreibung in JAGS ist fast identisch mit der mathematischen Beschreibung:

• “~” steht für “verteilt wie”
• dpois bzw. dnorm steht für Poisson-Verteilung bzw. Normalverteilung
• Vorsicht bei dnorm: Der erste Parameter ist wie erwartet der Mittelwert, zweiter Parameter ist jedoch nicht die Varianz, sondern die Präzision, also die inverse Varianz!
• statt einem Index wie \(y_t\) benutzen wir eckige Klammer \(y[t]\)
• Wir haben \(T\) Daten (Zeitpunkte), also müssen wir \(y_1, y_2, \ldots, y_T\) definieren; dazu können wir for-Schleifen wie in R benutzen
• log(lambda[t]) berechnen wir aus alpha und beta, wie in R benutzen wir “<-” als Zuweisung
• # bezeichnet einen Kommentar (wie in R)

Damit ist das Modell schon fertig beschrieben.

rjags arbeitet eigentlich mit Text-Dateien. Daher definieren wir hier die Variable model als Text (Anführungszeichen am Anfang und am Ende).

Nebenbei bemerkt: Modell meint hier sowohl das Datenmodell als auch die Priori. Wie schon bei der Regularisierungspriori gesehen, können wir auch mit den Prioris “modellieren”.

Daten

Neben dem Modell brauchen wir noch die Daten. Diese übergeben wir als Liste:

daten<-list(
  y=as.integer(y),
  e=as.integer(e),
  x=x,
  T=length(y)
)

• Die Daten müssen natürlich alles enthalten, was wir im Modell nicht genauer spezifiziert haben, hier also y, e, x, und T
• y und e wandeln wir in Integer um, also in ganze Zahlen

Kompilieren

Nun können wir das “Modell kompilieren”, sprich aus unserer Modellbeschreibung ein ausführbares JAGS-Modell machen:

jagsmodel <- jags.model(file = textConnection(model), data = daten)

## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 22
##    Unobserved stochastic nodes: 2
##    Total graph size: 159
## 
## Initializing model

jags.model() erwartet dabei eigentlich eine Datei als erstes Argument, wir übergeben hier den Text mittels textConnection().

Ziehen

Nun ziehen wir 1000 Mal aus dem Modell:

ziehungen <- coda.samples(jagsmodel, variable.names = c("alpha","beta"), n.iter=1000)

Mit dem plot-Befehl können wir uns Trace Plots und geschätzte Dichte der Posteriori für jeden Parameter erzeugen:

plot(ziehungen)

• Wie man an den Trace Plots sieht, sind die Ziehungen noch sehr instabil.
• Im Trace Plot sieht man hier zusätzlich einen laufenden Mittelwert (also Iteration gegen Mittelwert der Ziehungen bis zu dieser Iteration)
• Intern wird hier das coda-Paket benutzt, welches Ziehungen aus MCMC sinnvoll speichert und anzeigt.
• Die Iterationen beginnen hier bei 1001 - jags.model() hat nach der Kompilierung schon Mal 1000 Iterationen als burn-in gemacht.

Mehr Ziehungen

Ziehen wir also noch 20000 Mal:

ziehungen <- coda.samples(jagsmodel, variable.names = c("alpha","beta"), n.iter=20000)
plot(ziehungen)

Die Trace Plots sehen soweit gut aus (Nebenbemerkung: diese visuelle Inspektion ist natürlich subjektiv. Um das zu überprüfen, gibt es auch geeingete Maßzahlen).

Sehen wir uns also die Punkt- und Intervallschätzungen an:

summary(ziehungen)

## 
## Iterations = 2001:22000
## Thinning interval = 1 
## Number of chains = 1 
## Sample size per chain = 20000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##           Mean       SD  Naive SE Time-series SE
## alpha -0.04951 0.005831 4.123e-05      0.0002927
## beta   0.02536 0.002894 2.046e-05      0.0001425
## 
## 2. Quantiles for each variable:
## 
##           2.5%      25%      50%      75%    97.5%
## alpha -0.06089 -0.05351 -0.04950 -0.04560 -0.03770
## beta   0.01951  0.02342  0.02537  0.02733  0.03103

coda liefert uns hier:

1. Mittelwert (Mean) der Ziehungen (als Schätzer für den Posteriori-Erwartungswert) und Standardabweichung der Posterior (SD) für jeden Parameter (sowie geschätzter Fehler des Mittelwerts ohne (Naive SE) und mit (Time-series SE) Berücksichtigung der Autokorrelation der MCMC-Kette. Diese sollten möglichst klein sein)
2. (Posteriori-)Median und verschiedene Quantile (diese können natürlich auch gewählt werden)

Die Ergebnisse sind natürlich ganz ähnlich zu denen aus, die wir bei der Implementation von Hand hatten.

Ergänzen wir jetzt das Modell um die Überdispersion.

Einzige Änderungen am Modell:

• In der Gleichung für log(lambda) ergänzen wir + epsilon[t]
• Wir brauchen Prioris für epsilon[t] in der for-Schleife

model.overdisp <- "model{
    for (t in 1:T) {
        y[t] ~ dpois(lambda[t]*e[t])
        log(lambda[t]) <- alpha + beta*x[t] + epsilon[t]
        epsilon[t] ~ dnorm(0, 1000)
    }
    ## Prioris
    alpha ~ dnorm(0, 0.001)
    beta ~ dnorm(0, 0.001)
}"

Daten

Die Daten können wir von oben übernehmen.

Startwerte

• Unsere MCMC-Algorithmus braucht eigentlich Startwerte. Diese hatten wir zuvor nicht spezifiziert; JAGS zieht sie aus der Priori!
• Wir können sie aber auch (als Liste) vorgeben, was in der Regel Rechenzeit spart:

startwerte <- list(
  "alpha" = 0,
  "beta" = 0,
  "epsilon" = rep(0,T)
)

jagsmodel.overdisp <- jags.model(file = textConnection(model.overdisp), 
                                 data = daten, inits = startwerte)

## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 22
##    Unobserved stochastic nodes: 24
##    Total graph size: 182
## 
## Initializing model

Ziehungen

Wir ziehen erst 50000 Mal und verwerfen die Ziehungen als burn-in:

update(jagsmodel.overdisp, n.iter = 50000)

Dann die Ziehungen, bei denen wir mitprotokollieren. variable.names gibt an, welche Parameter wir abspeichern (später werden uns vielleicht nicht mehr alle Parameter interessieren):

ziehungen <- coda.samples(jagsmodel.overdisp, variable.names = c("alpha","beta", "epsilon"), n.iter=50000)
summary(ziehungen)

## 
## Iterations = 51001:101000
## Thinning interval = 1 
## Number of chains = 1 
## Sample size per chain = 50000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##                   Mean       SD  Naive SE Time-series SE
## alpha       -4.911e-02 0.028238 1.263e-04      0.0041154
## beta         2.550e-02 0.014047 6.282e-05      0.0021200
## epsilon[1]  -3.085e-02 0.015820 7.075e-05      0.0015526
## epsilon[2]  -1.727e-02 0.015037 6.725e-05      0.0013803
## epsilon[3]  -6.731e-03 0.013883 6.208e-05      0.0011359
## epsilon[4]   1.959e-02 0.012704 5.682e-05      0.0008832
## epsilon[5]  -5.218e-03 0.011968 5.352e-05      0.0007324
## epsilon[6]   1.271e-02 0.011101 4.965e-05      0.0005695
## epsilon[7]   2.788e-03 0.010631 4.754e-05      0.0004688
## epsilon[8]   4.467e-03 0.009741 4.356e-05      0.0002990
## epsilon[9]   1.211e-02 0.009485 4.242e-05      0.0002530
## epsilon[10]  1.377e-03 0.009475 4.237e-05      0.0002376
## epsilon[11]  3.035e-02 0.009525 4.260e-05      0.0002427
## epsilon[12]  9.060e-03 0.009805 4.385e-05      0.0002623
## epsilon[13]  2.031e-02 0.010122 4.527e-05      0.0003028
## epsilon[14]  3.369e-02 0.010314 4.612e-05      0.0003325
## epsilon[15]  2.243e-02 0.010831 4.844e-05      0.0003956
## epsilon[16]  1.209e-02 0.011098 4.963e-05      0.0004291
## epsilon[17] -4.698e-05 0.011135 4.980e-05      0.0004301
## epsilon[18] -1.159e-02 0.011268 5.039e-05      0.0004534
## epsilon[19] -2.649e-02 0.011392 5.095e-05      0.0004496
## epsilon[20] -3.180e-02 0.011379 5.089e-05      0.0004477
## epsilon[21] -2.573e-02 0.011240 5.027e-05      0.0004499
## epsilon[22] -3.344e-02 0.011119 4.973e-05      0.0004341
## 
## 2. Quantiles for each variable:
## 
##                   2.5%       25%        50%       75%      97.5%
## alpha       -0.1073074 -0.067911 -4.978e-02 -0.030307  0.0054134
## beta        -0.0017622  0.016124  2.573e-02  0.034917  0.0537035
## epsilon[1]  -0.0617562 -0.041574 -3.073e-02 -0.020275  0.0003858
## epsilon[2]  -0.0467995 -0.027426 -1.721e-02 -0.007253  0.0127450
## epsilon[3]  -0.0342224 -0.016041 -6.602e-03  0.002608  0.0205943
## epsilon[4]  -0.0051737  0.010953  1.965e-02  0.028125  0.0447592
## epsilon[5]  -0.0287219 -0.013264 -5.241e-03  0.002803  0.0186436
## epsilon[6]  -0.0088388  0.005175  1.264e-02  0.020096  0.0347859
## epsilon[7]  -0.0179267 -0.004396  2.814e-03  0.009837  0.0238559
## epsilon[8]  -0.0145083 -0.002124  4.391e-03  0.011005  0.0236518
## epsilon[9]  -0.0063129  0.005676  1.210e-02  0.018498  0.0307118
## epsilon[10] -0.0169066 -0.005088  1.332e-03  0.007758  0.0200136
## epsilon[11]  0.0117830  0.023951  3.029e-02  0.036751  0.0490554
## epsilon[12] -0.0101215  0.002443  9.013e-03  0.015709  0.0282944
## epsilon[13]  0.0005794  0.013496  2.030e-02  0.027128  0.0402292
## epsilon[14]  0.0136403  0.026669  3.367e-02  0.040673  0.0539005
## epsilon[15]  0.0013427  0.015094  2.236e-02  0.029623  0.0439424
## epsilon[16] -0.0094754  0.004568  1.206e-02  0.019577  0.0338201
## epsilon[17] -0.0216074 -0.007615 -9.971e-05  0.007480  0.0217964
## epsilon[18] -0.0334334 -0.019189 -1.162e-02 -0.004030  0.0104875
## epsilon[19] -0.0483239 -0.034274 -2.654e-02 -0.018816 -0.0040228
## epsilon[20] -0.0537859 -0.039467 -3.194e-02 -0.024155 -0.0091982
## epsilon[21] -0.0477025 -0.033317 -2.574e-02 -0.018148 -0.0037717
## epsilon[22] -0.0551418 -0.040887 -3.349e-02 -0.026023 -0.0116127

Trace Plots und Dichteschätzung von alpha

Wir plotten hier nicht alle Trace Plots, sondern beschränken uns auf alpha

plot(ziehungen[[1]][,1])

Trace Plots und Dichteschätzung von beta

… und beta:

plot(ziehungen[[1]][,2])

Sehen wir uns zuletzt noch die geschätzten Überdispersionsparameter an:

epsilon.sample <- ziehungen[[1]][,3:24]
epsilon.median <- apply(epsilon.sample, 2, median)
epsilon.q <- apply(epsilon.sample, 2, quantile, c(.025,.975))
plot(epsilon.median, ylim=range(epsilon.q), ylab=expression(epsilon))
lines(epsilon.q[1,])
lines(epsilon.q[2,])

Wir sehen wieder eine Art Trend, wenn auch etwas anders als zuvor geschätzt.

• Bisher hatten wir uns die Priori-Varianz der Überdispersion \(\epsilon\) subjektiv vorgegeben.
• Können wir diese aus den Daten schätzen?

Hierarchisches Modell

• Nehmen wir also als Priori für \(\epsilon_t\) an

\[ \epsilon_t \sim N(0,1/\tau) \]

• \(\tau\) ist die unbekannte inverse Varianz (Präzision; wir könnten auch mit der Varianz selbst arbeiten, Bayesianer nehmen aber lieber Präzision).
• “Aus den Daten schätzen” heißt: Nach der Beobachtung haben wir Posteriori-Information über \(\tau\).
• Wir brauchen also auch Priori-Information über \(\tau\).
• Wir nehmen uns ein Beispiel an der Priori für die Varianz der Normalverteilung, also

\[ \tau \sim Ga(a,b) \]

Die Präzision \(\tau\) ist damit Gamma-verteilt, was identisch ist zu die Varianz \(1/\tau\) ist Invers-Gamma-verteilt. \(a\) und \(b\) sind Hyperparameter, die wir noch wählen müssen.

Änderungen am Modell

• Wir ergänzen tau in der Priori der epsilon[t]
• Wir ergänzen die Priori von tau
• Wir wählen hier \(a=1\) und \(b=0.001\), der Priori-Erwartungswert ist damit 1000

model.overdisp2 <- "model{
    for (t in 1:T) {
        y[t] ~ dpois(lambda[t]*e[t])
        log(lambda[t]) <- alpha + beta*x[t] + epsilon[t]
        epsilon[t] ~ dnorm(0, tau)
    }
    ## Prioris
    alpha ~ dnorm(0, 0.001)
    beta ~ dnorm(0, 0.001)
    tau ~ dgamma(1, 0.001) 
}"

Startwerte und Kompilieren

• Spezifieren wir noch Startwerte:

startwerte <- list(
  "alpha" = 0,
  "beta" = 0,
  "tau" = 1000,
  "epsilon" = rep(0,T)
)

• und kompilieren das Modell:

jagsmodel.overdisp2 <- jags.model(file = textConnection(model.overdisp2), 
                                 data = daten, inits = startwerte)

## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 22
##    Unobserved stochastic nodes: 25
##    Total graph size: 183
## 
## Initializing model

Ziehungen

• 50000 Ziehungen burn-in (das ist übrigens sehr großzügig)

update(jagsmodel.overdisp2, n.iter = 50000)

• und schliesslich die endgültigen Ziehungen

ziehungen <- coda.samples(jagsmodel.overdisp2, variable.names = c("alpha","beta", "tau", "epsilon"), n.iter=50000)
summary(ziehungen)

## 
## Iterations = 51001:101000
## Thinning interval = 1 
## Number of chains = 1 
## Sample size per chain = 50000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##                   Mean        SD  Naive SE Time-series SE
## alpha       -5.106e-02 2.275e-02 1.018e-04      0.0026423
## beta         2.630e-02 1.142e-02 5.108e-05      0.0013752
## epsilon[1]  -2.842e-02 1.291e-02 5.776e-05      0.0010088
## epsilon[2]  -1.543e-02 1.224e-02 5.474e-05      0.0008847
## epsilon[3]  -5.630e-03 1.148e-02 5.132e-05      0.0007790
## epsilon[4]   1.958e-02 1.053e-02 4.707e-05      0.0006036
## epsilon[5]  -4.298e-03 9.902e-03 4.428e-05      0.0004854
## epsilon[6]   1.285e-02 9.383e-03 4.196e-05      0.0003911
## epsilon[7]   3.279e-03 9.045e-03 4.045e-05      0.0002950
## epsilon[8]   4.703e-03 8.409e-03 3.761e-05      0.0001791
## epsilon[9]   1.200e-02 8.194e-03 3.665e-05      0.0001530
## epsilon[10]  1.678e-03 8.250e-03 3.689e-05      0.0001485
## epsilon[11]  2.932e-02 8.269e-03 3.698e-05      0.0001485
## epsilon[12]  8.845e-03 8.444e-03 3.776e-05      0.0001689
## epsilon[13]  1.959e-02 8.754e-03 3.915e-05      0.0002015
## epsilon[14]  3.245e-02 8.886e-03 3.974e-05      0.0002391
## epsilon[15]  2.157e-02 9.302e-03 4.160e-05      0.0002924
## epsilon[16]  1.162e-02 9.475e-03 4.237e-05      0.0003186
## epsilon[17] -4.777e-05 9.461e-03 4.231e-05      0.0003249
## epsilon[18] -1.106e-02 9.538e-03 4.266e-05      0.0003370
## epsilon[19] -2.542e-02 9.619e-03 4.302e-05      0.0003323
## epsilon[20] -3.051e-02 9.681e-03 4.329e-05      0.0003426
## epsilon[21] -2.471e-02 9.478e-03 4.239e-05      0.0003368
## epsilon[22] -3.203e-02 9.449e-03 4.226e-05      0.0003127
## tau          2.020e+03 6.642e+02 2.970e+00      8.7742640
## 
## 2. Quantiles for each variable:
## 
##                   2.5%        25%        50%        75%      97.5%
## alpha        -0.094161 -6.650e-02 -5.087e-02 -3.714e-02 -4.396e-03
## beta          0.003232  1.938e-02  2.601e-02  3.429e-02  4.835e-02
## epsilon[1]   -0.055269 -3.668e-02 -2.813e-02 -1.972e-02 -4.088e-03
## epsilon[2]   -0.040336 -2.336e-02 -1.519e-02 -7.165e-03  8.159e-03
## epsilon[3]   -0.028890 -1.311e-02 -5.549e-03  2.157e-03  1.640e-02
## epsilon[4]   -0.001568  1.264e-02  1.962e-02  2.665e-02  3.992e-02
## epsilon[5]   -0.024382 -1.079e-02 -4.197e-03  2.378e-03  1.486e-02
## epsilon[6]   -0.005831  6.588e-03  1.293e-02  1.919e-02  3.105e-02
## epsilon[7]   -0.014784 -2.717e-03  3.384e-03  9.342e-03  2.085e-02
## epsilon[8]   -0.011973 -8.977e-04  4.722e-03  1.036e-02  2.108e-02
## epsilon[9]   -0.003990  6.483e-03  1.199e-02  1.753e-02  2.810e-02
## epsilon[10]  -0.014426 -3.896e-03  1.662e-03  7.282e-03  1.777e-02
## epsilon[11]   0.013267  2.377e-02  2.926e-02  3.485e-02  4.571e-02
## epsilon[12]  -0.007742  3.215e-03  8.791e-03  1.451e-02  2.548e-02
## epsilon[13]   0.002657  1.369e-02  1.953e-02  2.542e-02  3.700e-02
## epsilon[14]   0.015349  2.646e-02  3.234e-02  3.833e-02  5.020e-02
## epsilon[15]   0.003222  1.536e-02  2.162e-02  2.774e-02  3.994e-02
## epsilon[16]  -0.007127  5.325e-03  1.160e-02  1.790e-02  3.016e-02
## epsilon[17]  -0.018625 -6.341e-03  2.134e-05  6.160e-03  1.849e-02
## epsilon[18]  -0.029806 -1.731e-02 -1.104e-02 -4.676e-03  7.568e-03
## epsilon[19]  -0.044674 -3.161e-02 -2.539e-02 -1.900e-02 -6.606e-03
## epsilon[20]  -0.049785 -3.688e-02 -3.038e-02 -2.408e-02 -1.159e-02
## epsilon[21]  -0.043508 -3.093e-02 -2.461e-02 -1.838e-02 -6.285e-03
## epsilon[22]  -0.050942 -3.820e-02 -3.194e-02 -2.578e-02 -1.367e-02
## tau         961.392461  1.546e+03  1.941e+03  2.407e+03  3.549e+03

• Insbesondere erhalten wir hier einen geschätzten Wert von \(\tau\) von 2046 (Posteriori-Erwartungswert) bzw. 1971 (Posteriori-Median)
• Posteriori-Erwartungswert und Posteriori-Median unterscheiden sich hier mehr als vergleichsweise bei alpha und beta: Die Posteriori-Dichte ist nämlich nicht symmetrisch.

Trace Plots und Dichteschätzung von alpha, beta und tau

plot(ziehungen[[1]][,1])

plot(ziehungen[[1]][,2])

plot(ziehungen[[1]][,25])

Geschätzten Überdispersionsparameter

epsilon.sample <- ziehungen[[1]][,3:24]
epsilon.median <- apply(epsilon.sample, 2, median)
epsilon.q <- apply(epsilon.sample, 2, quantile, c(.025,.975))
plot(epsilon.median, ylim=range(epsilon.q), ylab=expression(epsilon))
lines(epsilon.q[1,])
lines(epsilon.q[2,])

Wir sehen wieder eine Art Trend, wenn auch etwas anders als zuvor geschätzt.

• JAGS ermöglicht es uns, Bayesianische Modelle einfach zu schätzen.
• Im letzten Beispiel hatten wir eine Priori-Verteilung auf einen Parameter in einer Priori (\(\tau\)). Man nennt diese Hyper-Priori
• Da die Hyper-Priori auf einer anderen Ebene des Modells arbeitet (nicht direkt auf den Parametern der Datendichte), nennt man diese Art von Modell auch Hierarchisches Modell (dazu später mehr)

• Weiter zum Quiz
• Zurück zur Übersicht