Skip to Tutorial Content

Drivers-Beispiel

Daten

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

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

\[ \gamma_i|\gamma_{i-1},\tau_c \sim N(\gamma_{i-1},\tau_c^{-1}); j=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.\).
  • Daraus ergibt sich für die gemeinsame Verteilung der \(\gamma\):

\[ p(\gamma|\tau_c)\propto \tau_c^{T-1} \exp\left(-\frac{\tau_c}{2}\sum_{i=2}^T(\gamma_i-\gamma_{i-1})^2\right) \]

Dabei handelt es sich um eine uneigentliche Normalverteilung. Das heißt: Die Dichte hat die Form einer Normalverteilung, aber es gilt:

\[ \int p(\gamma|\tau_c)d\gamma = +\infty \]

Warum ist das so?

  • Für unsere \(T\)-dimensionale Zufallsvariable haben wir \(T-1\) “richtige” Normalverteilungen als Randverteilungen definiert.

  • Für \(\gamma_1)\) ab haben wir \(p(\gamma_1)\propto const.\), also eine uneigentliche (Gleich-)Verteilung.

  • Es fehlt also für die \(T\)-dimensionale Verteilung eine Dimension.

  • Wir behandeln diese uneigentliche Verteilung aber trotzdem wie eine übliche mehrdimensionale Normalverteilung.

  • Wir wissen aus den einfachen Fällen, dass wir uneigentliche Prioris benutzen dürfen, nur die Posteriori mur eine eigentliche Verteilung sein.

Matrixschreibweise

Wir gehen für die mehrdimensionale Verteilung zur Matrixschreibweise über.

Sei \(\gamma=(\gamma_1,\gamma_2,\ldots,\gamma_T)\) und \(\Delta(\gamma)=((\gamma_2-\gamma_1),(\gamma_3-\gamma_2),\ldots,(\gamma_T-\gamma_{T-1}))\) der Vektor der Differenzen. Dann kann man eine Differenzenmatrix definieren, die aus \(\gamma\) \(\Delta(\gamma)\) macht:

\[ D = \left(\begin{array}{rrrrrr} -1 & 1 & 0 & & \cdots & 0\\ 0 & -1 & 1 & 0 & \cdots & 0\\ \vdots & & \ddots& \ddots & & \vdots\\ 0 & & \cdots & 0 & -1 & 1 \end{array}\right) \]

Also ist \(D\cdot \gamma=\Delta(\gamma)\).

Nun quadrieren wir \(\Delta(\gamma)\) und summieren auf:

\[ \begin{aligned} \Delta(\gamma)'\Delta(\gamma)&=\sum_{i=2}^T(\gamma_i-\gamma_{i-1})^2\\ \Delta(\gamma)'\Delta(\gamma)&=(D\cdot \gamma)'(D\cdot \gamma)=\gamma'(D'D)\gamma\\ \end{aligned} \]

mit

\[ D'D= Q_c = \left(\begin{array}{rrrrrr} 1 & -1 & & & & \\ -1 & 2 & -1 & & & \\ & -1 & 2 & -1 & & \\ & & \ddots & \ddots & \ddots & \\ & & & -1 & 2 & -1\\ & & & & -1 & 1 \end{array}\right) \]

Insgesamt ist also

\[ p(\gamma|\tau_c)\propto \tau_c^{T-1} \exp\left(-\frac{\tau_c}{2}\gamma'Q_c\gamma\right) \]

Das sieht nach der Dichte einer \(T\)-dimensionalen Normalverteilung aus:

\[ \gamma \sim N(0,\tau_cQ_c^{-1}) \]

\(\tau_cQ_c^{-1}\) ist dabei die inverse Kovarianzmatrix (Präzisionsmatrix) der Normalverteilung.

  • Leider ist \(Q_c\) gar nicht invertierbar (erkennt man z.B. daran, dass alle Zeilensummen gleich Null sind)
  • Das hängt damit zusammen, dass wir nur eine uneigentliche Normalverteilung haben. Diese hat keine Varianz und auch keinen Mittelwert.
  • Die Präzisionsmatrix ist aber definiert und mit dieser können wir weiterarbeiten

Das Matrix-R-Paket erlaubt es, die Präzisionsmatrix als spärliche (sparse) Matrix zu speichern, das heißt, nur Einträge ungleich Null werden abgespeichert.

library(Matrix)
Qc <- Matrix::sparseMatrix(i = c(1:T,1:(T-1)),
        j=c(1:T,2:T), 
        x=c(1,rep(2,T-2),1,rep(-1,T-1)),
        symmetric = TRUE)
image(Qc)

Man erkennt gut die “Bandstruktur” der Matrix: Nur die Hautdiagonale und die erste Nebendiagonale ist besetzt.

Priorimodell - Saisontrend

  • Nun wollen wir den Saisontrend modellieren.
  • Möglich wäre z.B. parametrische Sinus-Kurve, das wäre vermutlich zu unflexibel
  • Stattdessen benutzten wir wieder Vorwissen, wie der Saisontrend aussieht:

\[ \begin{aligned} p(\delta|\tau_d)&\propto (\tau_d^{T-12+1})/2\cdot\\ &\cdot\exp\left(-\frac{\tau_d}{2}\sum_{i=0}^{T-12}(\delta_{i+1}+\delta_{i+2}+\ldots+\delta_{i+12})^2\right) \end{aligned} \]

  • Es lässt sich zeigen:

\[ \delta|\tau_d \sim N_T\left(0,(\tau_d Q_d)^{-1}\right) \]

mit

\[ Q_d=\left(\begin{array}{cccccccccccccc} 1 & 1 & 1 & & \cdots && 1 & 1 & 0 & &\cdots & & 0 \\ 1 & 2 & 2 & & \cdots & & 2 & 2 & 1 &0 & \cdots & & 0\\ 1 & 2 & 3 & & \cdots & & 3 & 3 & 2 & 1 & \ddots & \\ \vdots & \vdots & \vdots & & \ddots \\ 1 & 2 & 3 & \cdots & 11 & 12 &\cdots & 12& 11 & 10 & \cdots & & 0\\ 0 & 1 & 2 & 3 & \cdots & 11 & 12 &\cdots & 12& 11 & 10 & \cdots &0\\ 0 & \ddots& \ddots & \ddots & \ddots & \ddots & \ddots & \ddots& \ddots & \ddots & \ddots & \\ \end{array}\right) \]

Qd <- Matrix::sparseMatrix(i=T,j=T,x=0,
        symmetric = TRUE)
EinsM<-Matrix::sparseMatrix(i=rep(1:12,12),
        j=rep(1:12,each=12), x=1)
for (i in 1:(T-12))
  Qd[i+(0:11),i+(0:11)]<-Qd[i+(0:11),i+(0:11)]+EinsM
image(Qd)

Auch hier erkennt man die “Bandstruktur” der Matrix, mit einem breiteren Band: Die Hautdiagonale und die ersten elf Nebendiagonale sind besetzt, der Rest der Einträge ist 0.

Priori-Modell - Intercept und Kovariableneffekt

  • Keine Priori-Information

\[p(\alpha)\propto const. \Leftrightarrow ''\alpha\sim N\left(0,0^{-1}\right)'' \]

  • Analog \(p(\beta)\propto const.\)
  • Also: Alle Effekte a priori normalverteilt, enthalten unbekannte Präzisionsparameter \(\tau_c\) bzw. \(\tau_d\)
  • Gemeinsame Priori von \(\theta=(\alpha,\beta,\gamma,\delta)\) hat Form

\[ \begin{aligned} \theta|\tau &\sim N_{2+2T}(0,Q(\tau)^{-1})\\ Q(\tau)&=\left(\begin{array}{cccc} 0&0&0&0\\ 0&0&0&0\\ 0&0&\tau_cQ_c &0\\ 0&0&0&\tau_dQ_d\\ \end{array}\right) \end{aligned} \]

Hyperpriors

  • Die Präzisionsparameter in den Prioris steuern, wie stark die Regularisierung des Zeit- bzw. Saisontrends ist
  • Diese Parameter der Prioris, Hyperparameter, wollen wir nicht vorgeben, sondern mitschätzen
  • Jeder Parameter, der geschätzt werden soll, braucht eine Priori, hier Hyperpriori genannt

\[ \begin{aligned} \tau_c\sim Ga(a_c,b_c)\\ \tau_d\sim Ga(a_d,b_d) \end{aligned} \]

  • Außerdem brauchen wir noch eine Priori für \(\sigma^2\sim IG(a_s,b_s)\)

Posteriori

Sei \(y^*=\sqrt{y}\)

\[ \begin{aligned} p(\theta,\tau|y^*) &\propto f(y^*|\theta) p(\theta,\tau)\\ &\propto f(y^*|\theta) p(\theta|\tau) p(\tau) \\ &\propto \prod_{i=1}^T \sigma^{-1} \exp\left(-\frac{1}{2\sigma^2}(y_i^*-\alpha-\beta x_i-\gamma_i+\delta_i )^2\right) \cdot\\ &\cdot \tau_c^{(T-1)/2} \exp\left(-\frac{\tau_c}{2}\gamma^TQ_c\gamma\right) \tau_d^{(T-m+1)/2} \exp\left(-\frac{\tau_d}{2}\delta^TQ_d\delta\right)\\ &\cdot \tau_c^{a_c-1}\exp(-\tau_c b_c)\cdot \tau_d^{a_d-1}\exp(-\tau_d b_d)\\ &\cdot (\sigma^2)^{-a_s-1}\exp(-b_s/\sigma^2) \end{aligned} \]

Full conditional

Zeitlicher Trend

\[ \begin{aligned} p(\gamma|\cdot) &\propto \exp\left(-\frac{1}{2\sigma^2}\sum_i(y_i^*-\gamma_i-\delta_i-\beta x_i-\alpha )^2\right) \cdot \exp\left(-\frac{\tau_c}{2}\gamma^TQ_c\gamma\right)\\ &\propto\exp\left(-\frac{1}{2\sigma^2}\sum_i(\gamma_i-\epsilon^{c}_i)^2-\frac{\tau_c}{2}\gamma^TQ_c\gamma\right)\\ &\propto \exp\left(-\frac{1}{2}\gamma^T\frac{T}{\sigma^2}I\gamma+\frac{1}{\sigma^2}\gamma^T\epsilon_c^*-\frac{1}{2}\gamma^T\tau_cQ_c\gamma\right)\\ &\propto \exp\left(-\frac{1}{2}\gamma^T\left(\frac{T}{\sigma^2}I+\tau_cQ_c\right)\gamma+\frac{1}{\sigma^2}\gamma^T\epsilon_c^*\right) \end{aligned} \]

Kanonische Form der Normalverteilung mit Präzisionsmatrix \(Q_c^*=\left(\frac{1}{\sigma^2}I+\tau_cQ_c\right)\) und Erwartungswert \(Q_c^{*-1}m_c\) mit \(m_{c,i}=\frac{1}{\sigma^2}\sum_{i}(y_i^*-\alpha-\beta x_i-\delta_i)\).

Saison-Effekt

Analog: \[ \begin{aligned} \delta|\cdot &\sim N_{T}((Q_d^*)^{-1}m_d,(Q_d^*)^{-1})\\ Q_d^* &=\left(\frac{1}{\sigma^2}I+\tau_dQ_d\right)\\ m_{d,k}&=\frac{1}{\sigma^2}\sum_{i}(y_i^*-\alpha-\beta x_i-\gamma_i) \end{aligned} \]

Kovariableneffekt

\[ \begin{aligned} p(\beta|\cdot) &\propto \exp\left(-\frac{1}{2\sigma^2}\sum_i(y_i^*-\gamma-\delta-\alpha-\beta x_i)^2\right)\\ &\propto \exp\left(-\frac{1}{2}\beta^2\frac{1}{\sigma^2}\sum_i x_i^2 + \beta \frac{1}{\sigma^2}\sum_i(x_iy_i^*-x_i\alpha-x_i\gamma_i-x_i\delta_i)\right) \end{aligned}\]

Also \(\beta|\cdot \sim N(q_b^{-1}m_b,q_b^{-1})\) mit \(q_b=\frac{\sum x_i^2}{\sigma^2}\) und \(m_b=\frac{1}{\sigma 2}\sum_i(x_iy_i^*-x_i\alpha-x_i\gamma_j-x_i\delta_k)\)

  • Analog: \(\alpha|\cdot \sim N(q_a^{-1}m_a,q_a^{-1})\) mit \(q_a=\frac{T}{\sigma^2}\) und \(m_b=\frac{1}{\sigma 2}\sum_i(y_i^*-x_i\beta-\gamma_j-\delta_k)\)

Präzisionsparameter

Allgemein: \(p(\tau|\cdot) \propto p(\theta|\tau)p(\tau)\), also ist die Full Conditional unabhängig von den Daten

Konkret:

\[ \begin{aligned} p(\tau_c|\cdot) &\propto \tau_c^{(T-1)/2} \exp{\left(\frac{\tau_c}{2}\gamma^TQ_c\gamma\right)} \tau_c^{a_c-1}\exp(-\tau_c b_c)\\ &\propto \tau_c^{c+(T-1)/2-1}\exp\left(-\tau_c(b_c+\gamma^T Q_c\gamma/2) \right) \end{aligned} \]

Also

\[ \begin{aligned} \tau_c|\cdot &\sim Ga\left(a_c+(T-1)/2, b_c+\gamma^T Q_c\gamma/2\right)\\ \tau_d|\cdot &\sim Ga\left(a_d+(T-m+1)/2, b_d+\delta^T Q_d\delta/2\right)\\ \sigma_2|\cdot &\sim IG\left(a_s+T/2, b_s+\sum(y_i^*-\alpha-x_i\beta+\gamma_j-\delta_k)^2\right) \end{aligned}\]

Implementation

Vorarbeiten

y <- sqrt(Drivers$y)
belt <- Drivers$belt

a.c <- a.d <- 1
b.c <- 0.0005
b.d <- 0.1
a.s <- 1
b.s <- 1

alpha <- mean(y)
beta <- 0
gamma <- rep(0,T)
delta <- rep(0,T)
tau.c <- a.c/b.c
tau.d <- a.d/b.d
sigma2 <- b.s/a.s
sumx2 <- sum(belt)
burnin=2000
nr.it=10000
I = burnin+nr.it

alpha.save<-beta.save<-sigma2.save<-rep(0,nr.it)
gamma.save<-array(0,c(T,nr.it))
delta.save<-array(0,c(T,nr.it))
tau.save<-array(0,c(2,nr.it))
iter=0

Gibbs-Sampler

while (iter<I)
{
  iter=iter+1
  gamma<-gamma-mean(gamma);  delta<-delta-mean(delta)
  m <- sum(y-gamma-delta-belt*gamma)/sigma2; Q <- T/sigma2
  alpha <- rnorm(1, m/Q, 1/Q)
  m <- sum(belt*(y-gamma-delta-alpha))/sigma2; Q <- sumx2/sigma2
  beta <- rnorm(1, m/Q, 1/Q)
  m <- (y-delta-beta*belt-alpha)/sigma2; Q <- tau.c*Qc + Matrix::Diagonal(T)/sigma2
  gamma <- rwc::rnorm.Q(Q,as.matrix(m),canon=TRUE)[,1]
  m <- (y-gamma-beta*belt-alpha)/sigma2; Q <- tau.d*Qd + Matrix::Diagonal(T)/sigma2
  delta <- rwc::rnorm.Q(Q,as.matrix(m),canon=TRUE)[,1]
  tau.c <- rgamma(1, a.c + (T-1)/2, b.c + (t(gamma)%*%Qc%*%gamma/2)[1,1])
  tau.d <- rgamma(1, a.d + (T-12+1)/2, b.d + (t(delta)%*%Qd%*%delta/2)[1,1])
  sigma2 <- 1/rgamma(1, a.s + T/2, b.s + sum((y-alpha-gamma-delta-beta*belt)^2))
  
  if (iter>burnin){alpha.save[iter-burnin]<-alpha;  beta.save[iter-burnin]<-beta
    gamma.save[,iter-burnin]<-gamma; delta.save[,iter-burnin]<-delta
    tau.save[,iter-burnin]<-c(tau.c,tau.d); sigma2.save[iter-burnin]<-sigma2
  }
}

Ergebnisse

Intercept

summary(alpha.save)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   40.65   40.86   40.91   40.91   40.95   41.20
quantile(alpha.save, c(0.025,0.975))
##     2.5%    97.5% 
## 40.77421 41.03779

Kovariableneffekt

summary(beta.save)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  -5.181  -3.435  -3.091  -3.092  -2.757  -1.037
quantile(beta.save, c(0.025,0.975))
##      2.5%     97.5% 
## -4.068885 -2.114143

Prioriparameter

plot.ts(t(tau.save))

plot(density(tau.save[1,]))

plot(density(tau.save[2,]))

apply(tau.save, 1, summary)
##               [,1]      [,2]
## Min.      3.500446   4.47620
## 1st Qu.  10.509714  23.14163
## Median   15.143450  31.72772
## Mean     18.150302  34.37229
## 3rd Qu.  21.681792  42.83820
## Max.    175.769277 142.26650
apply(tau.save, 1, quantile, c(0.025,0.975))
##           [,1]     [,2]
## 2.5%   5.43638 12.51821
## 97.5% 47.79207 71.89869

Varianz

plot.ts(sigma2.save)

plot(density(sigma2.save))

summary(sigma2.save)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   2.928   4.440   4.890   4.964   5.400   8.848
quantile(sigma2.save, c(0.025,0.975))
##     2.5%    97.5% 
## 3.708657 6.620097

Zeiteffekt

gamma<-apply(gamma.save,1,quantile,c(.025,.5,.975))
plot(gamma[2,],type="l",ylim=range(gamma),
     xlab="Zeit in Monaten", ylab="", axes=FALSE)
axis(2, las=1)
axis(1,at=1+12*(0:16),labels = 1969:1985)
lines(gamma[1,],col=grey(.7))
lines(gamma[3,],col=grey(.7))
title(main="Zeiteffekt")

Saisoneffekt

delta<-apply(delta.save,1,quantile,c(.025,.5,.975))
plot(delta[2,],type="l",ylim=range(delta),
     xlab="Zeit in Monaten", ylab="", axes=FALSE)
axis(2, las=1)
axis(1,at=1+12*(0:16),labels = 1969:1985)
lines(delta[1,],col=grey(.7))
lines(delta[3,],col=grey(.7))
title(main="Saisoneffekt")

Glatte Zeitreihe

y.save<-gamma.save+delta.save
for (i in 1:T)
  y.save[i,]<-y.save[i,]+alpha.save+beta.save*belt[i]
y.save<-y.save^2

Weiter

Drivers-Beispiel