• Mehr Parameter führen automatisch zu besserer Anpassung des Modells an die Daten.

• Übliche Informationskriterien zur Modellwahl (z.B. AIC und BIC) messen daher die Anpassung des Modells an die Daten plus die Anzahl der Parameter.

• Bei hierarchischen Modellen ist die Anzahl der Parameter jedoch nicht aussagekräftig, da die Parameter schon a priori abhängig sind (siehe z.B. die Modellierung des Zeiteffekts)

• Bei hierarchischen Modellen ist es daher sinnvoll, die effektive Anzahl der Parameter zu schätzen.

• Oder: Wieviel Freiheit hat unser Modell wirklich?

Beispiel

• In unserem ersten hierarchischen Modell hatten wir für jeden Monat \(t=1,\ldots,T\) einen eigenen Zeiteffekt \(\gamma_t\).
• Damit hatten wir \(T\) \(\gamma-\)Parameter

Zusätzlich hatten wir als Priori-Annahme:

\[ \gamma_t|\gamma_{t-1},\tau_c \sim N(\gamma_{t-1},\tau_c^{-1}); t=2,\ldots,T \]

Für \(\tau_c\to\infty\), also \(\tau_c^{-1}\to0\) gilt

\[ \gamma_T=\gamma_{T-1}=\gamma_{T-2}=\ldots=\gamma_1 \]

• Alle Parameter sind dann gleich, also haben wir nur einen Parameter!
• Umgekehrt gilt für \(\tau_c\to0\), also \(\tau_c^{-1}\to\infty\), dass \(\gamma_1, \gamma_2,\ldots,\gamma_T\) (a prior) voneinander unbhängig sind. Dann haben wir wirklich \(T\) Parameter.
• In der Realität ist \(0<\tau_c<\infty\). Also ist die Anzahl der Parameter irgendwo zwischen 1 und \(T\)!
• Diese effektive Anzahl an Parametern wollen wir schätzen.

Der Ausdruck

\[ D(y) := -2 \Big( \log \big( p(y|\hat\theta)\big)-\log \big(p(y|\hat\theta_s)\big)\Big) \]

heißt Devianz.

Dabei sind

• \(\hat\theta\) die geschätzen Parameter eines Modells \(M\) und
• \(\hat\theta_s\) die geschätzten Parameter in einem saturierten Modell (Daten komplett angepasst).
• Dies entspricht minus zweimal dem Logarithmus des Likelihoodratios.
• Manchmal wird auch nur \(-2 \log \big( p(y|\hat\theta)\big)\) als Devianz bezeichnet.

Maß der Reduktion der Überraschung

Sei \(\theta^{TRUE}\) der wahre Parameter und \(\hat{\theta}\) ein Schätzer für \(\theta^{TRUE}\). Betrachte nun

\[ d(y,\theta^{TRUE},\hat{\theta})=-2\log\left(p(y|\theta^{TRUE})\right)+2\log\left(p(y|\hat{\theta})\right) \]

Dies kann interpretiert werden als Maß der Reduktion der Überraschung oder Unsicherheit aufgrund der Schätzung.

Vorschlag von Spiegelhalter et al.(2002): Benutze Posteriori-Erwartungswert von \(d(y,\theta^{TRUE},\hat{\theta})\) als Effektive Anzahl der Parameter im Bayes-Modell.

Effektive Anzahl von Parametern

\[ \begin{aligned} E(d(y,\theta^{TRUE},\hat{\theta})&=E(-2\log\left(p(y|\theta^{TRUE})\right))+E(2\log\left(p(y|\hat{\theta})\right))\\ &=: \widehat{D(\theta)}-D(\hat{\theta}) =: p_D \end{aligned} \]

\(p_D\) ist abhängig von

• den Daten \(y\)
• dem Datenmodell
• der Priori \(p(\theta)\)
• der Wahl des Punktschätzers \(\hat\theta\)

Das Deviance Information Criterion (DIC) ist definiert als

\[ DIC=\widehat{D(\theta)}+p_D \]

• Dies entspricht dem Bayesianischen Maß für die Anpassung des Fits plus den Komplexitätsterm \(p_D\).
• Es gilt: \(\widehat{D(\theta)}=E(D(\theta))=E(-2\log(p(y|\theta)))\) wird kleiner für bessere Anpassung.

• Bei substantiellem Konflikt zwischen Priori und Daten sowie bei nicht unimodalen Posterioris kann \(p_D\) negativ werden (dann nicht benutzen!).
• Falls kein hierarchische Modell vorliegt und alle Prioris komplett spezifiziert sind, gilt

\[ AIC=D(\hat{\theta}_{ML})+2\cdot p \]

Berechnung des DIC bei MCMC

• \(\widehat{D(\theta)}\): In jeder MCMC-Iteration \(k\) berechne \(D(\theta^{(k)}\) und schätze \(\hat{D(\theta)}\) über Mittelwert (Median)
• \({D(\hat\theta)}\): Benutze Punktschätzer (Mittelwert, Median) und setze diesen in \(D\) ein (plug-in Schätzer)

Poisson-Log-Normal-Modell

mu<-array(0,c(T,nr.it))
D <- rep(NA,nr.it)
for (i in 1:nr.it) # Pro Iteration
     { # Berechne Erwartungswert von y
        mu[,i]=alpha.save[i]+belt*beta.save[i]+
          gamma.save[,i]+delta.save[,i]
        # Devianz
        D[i] <- -2*sum(dnorm(y, mean=mu[,i], 
                    sd=sqrt(sigma2.save[i]),log=TRUE))
}

Wir bekommen also die Devianz (Anpassung) pro Iteration, insgesamt also die Posteriori der Devianz:

Dann berechnen wir die Devianz für den Punktschätzer \(D(\hat{\theta})\):

D.theta.hat <- -2*sum(dnorm(y,apply(mu,1,mean),
                  sqrt(mean(sigma2.save)),log=TRUE))

und vergleichen mit dem Punktschätzer der Devianz \(\hat{D}\)

D.hat <- mean(D)
pD <- D.hat - D.theta.hat

Alternativ kann man statt dem geschätzten Posteriori-Erwartungswert auch den Posteriori-Median nehmen:

D.theta.hat.med <- -2*sum(dnorm(y,apply(mu,1,median), sqrt(median(sigma2.save)),log=TRUE))
D.hat.med <- median(D)
pD.med <- D.hat.med - D.theta.hat.med

Die Werte unterscheiden sich in der Regel nur minimal.

Random Walk zweiter Ordnung

Als Vergleichsmodelle nehmen wir

• ein Modell komplett ohne Zeittrend und
• ein Modell nehmen wir einen Zeittrend mit Random Walk zweiter Ordnung als Priori (RW2):

\[ \gamma_i|\gamma_{i-1},\gamma_{i-2},\tau_c \sim N(2\gamma_{i-1}-\gamma_{i-2},\tau_c^{-1}); i=3,\ldots,T \]

• Das entspricht der Priori-Ananhme eines linearen Trends.
• Der Erwartungswert von \(\gamma_i\) ist \(\gamma_{i-1}+(\gamma_{i-1}-\gamma_{i-2})\).

Modellvergleich

• Hier hat also das Modell mit der RW2-Priori die beste Anpassung – die niedrigste Devianz –, geringfügig besser als das Modell mit RW1.
• Das Modell ohne Zeittrend hat die geringste Anzahl effektiver Parameter (\(p_D\)) – natürlich auch die geringste Anzahl real geschätzter Parameter, da kein Zeittrend geschätzt wird.
• Unter den Modell mit RW-Prioris hat RW1 geringfügig weniger effektive Parameter.
• Insgesamt ist damit das Modell mit RW2-Priori das beste, es hat den kleinsten DIC-Wert. Der Unterschied zum Modell mit RW1 ist jedoch nur gering. Das Modell ohne Zeittrend ist dagegen klar schlechter.
• Der DIC (ebenso wie Devianz und \(p_D\)) sind nur geschätzt, daher sollte man kleine Differenzen zwischen den Modellen nicht überinterpretieren. Konkret sind hier die Modelle mit RW1- und RW2-Priori praktisch gleich gut.