Beispiel MCMC bei Poisson-Lognormal

Als Beispiel sehen wir uns folgendes Modell an:

Uns liegen Daten über die Anzahl von Kaiserschnitten in 24 Krankenhäusern vor. Wir benutzen folgendes Modell:

\[ \begin{aligned} y_i&\sim Po(\lambda_i)\\ \log(\lambda_i)&\sim N(\mu,1)\\ \mu&\sim N(0,1000) \end{aligned} \]

Implementation

Unbekannte Parameter sind \(\mu\) und \(\tau\). Beide werden hier mittels Metropolis-Hastings-Schritten gezogen:

• Für \(\mu\) verwenden wir ein Random-Walk-Proposal: Auf den letzten Wert wird eine normalverteilte Zufallszahl addiert
• Für \(\tau\) verwenden wir eine grobe Approximation durch eine Normalverteilung.

poisson.lognormal.mcmc<-function(sumx,n,tau=1000,I=1000, do.tuning=FALSE, tuning=50, do.burnin=TRUE, burnin=100, rw.sd=.1){
mu<-rep(NA,I)
lambda<-rep(NA,I)
mu[1]<-1
lambda[1]<-sumx/n 

#Berechnung der log Full Conditional von mu
log.fc.mu <- function(mu,lambda,sumx,tau)
{
  return(mu*log(lambda)-0.5*mu^2-0.5*mu^2/tau)
}

log.fc.lambda <- function(lambda,mu,sumx,n){
  return((sumx+mu)*log(lambda)-n*lambda-0.5*log(lambda)^2)
}

akzeptanz <- 0

#Ziehen der Zufallszahlen mittel MH-Algorithmus
i=1
while (i<=I)
{
  i<-i+1
  
  # Ziehe mu mittels random walk proposal
  mustern <- rnorm(1,mu[i-1],rw.sd)
  logalpha<-log.fc.mu(mustern,lambda[i-1],sumx,tau)-log.fc.mu(mu[i-1],lambda[i-1],sumx,tau)
  alpha<-exp(logalpha)
  ifelse(runif(1)<alpha,
         {mu[i]<-mustern; akzeptanz=akzeptanz+1},
         mu[i]<-mu[i-1])
  
  # Ziehe lambda mittels ad hoc Normal-Approximation
  c<-sumx + mu[i]+1
  d<-n+1
  lambdastern<-rnorm(1,mean=c/d,sd=sqrt(c/d^2))
  logalpha <- log.fc.lambda(lambdastern,mu[i],sumx,n)-log.fc.lambda(lambda[i-1],mu[i],sumx,n)
  logalpha <- logalpha+dnorm(lambda[i-1],mean=c/d,sd=sqrt(c/d^2),log=TRUE)-dnorm(lambdastern,mean=c/d,sd=sqrt(c/d^2),log=TRUE)
  alpha<-exp(logalpha)
  ifelse (runif(1)<alpha,
          lambda[i]<-lambdastern,
          lambda[i]<-lambda[i-1])
  
  if (i==tuning)if(do.tuning)
  {
    # Tuning des RW-Proposals
    ak.rate = akzeptanz/tuning
    if (ak.rate>.5){
      rw.sd <-rw.sd * ak.rate * 2
      akzeptanz = 0
      i = 1
      lambda[1] <- lambda[tuning]
      mu[1] <- mu[tuning]
    }
    if (ak.rate<.3){
      rw.sd <-rw.sd * ak.rate
      akzeptanz = 0
      i = 1
      lambda[1] <- lambda[tuning]
      mu[1] <- mu[tuning]
    }
    message(paste("Akzeptanzrate Tuning:",round(ak.rate*100,1)))
    message(paste("Neuer Varianzparameter:",round(rw.sd,3),"\n"))
  }
  if(i==burnin)if (do.burnin)
  {
    i=1
    do.burnin=FALSE
    akzeptanz=0
  }
  
}
lambda<-lambda[-1]
mu<-mu[-1]
ak.rate = akzeptanz/i
message(paste("Akzeptanzrate gesamt:",round(ak.rate*100,1),"\n"))
return(coda::mcmc(cbind(mu,lambda)))
}

Durchführung

n <-dim(data)[1]
sumx <- sum(data$n)
mcmc.simple<-poisson.lognormal.mcmc(sumx,n)

## Akzeptanzrate gesamt: 96.7

Das coda-Paket stellt verschiedene Funktionen für die Analyse von MCMC-Ziehungen bereit. Im Beispiel wurden die Ziehungen als coda::mcmc-Objekt zurückgegeben.

summary

summary(mcmc.simple)

## 
## Iterations = 1:1000
## Thinning interval = 1 
## Number of chains = 1 
## Sample size per chain = 1000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##         Mean     SD Naive SE Time-series SE
## mu      2.72 0.6674  0.02111        0.24551
## lambda 10.52 0.6585  0.02082        0.04108
## 
## 2. Quantiles for each variable:
## 
##         2.5%    25%    50%    75%  97.5%
## mu     1.435  2.303  2.671  3.071  4.084
## lambda 9.315 10.064 10.477 10.985 11.872

plot

plot(mcmc.simple)

• An den trace plots ist erkennbar, dass die Ziehungen aus \(\lambda\) (vermutlich) “gut mischen”: Die Ziehungen decken relativ schnell den ganzen relevanten Wertebereicht der Posterioridichte ab.
• Dagegen mischt der Algorithmus von \(\mu\) relativ schlecht: Die Kette wandert langsam durch den Wertebereich. Das ist offensichtlich ineffizient.

Akzeptanzrate

Um das genauer zu beurteilen, kann man sich die Akzeptanzrate, gleich dem Anteil akzeptierter Vorschläge anschauen:

• Niedrige Akzeptanzrate: Kette bleibt oft im selben Zustand \(\to\) schlecht
• Zu hohe Akzeptanzrate: Kette bewegt sich (eventuell) nur langsam \(\to\) schlecht

Beispiel zu niedrige Akzeptranzrate

In obigem Beispiel ist die Akzeptanzrate (für \(\mu\)) fast 100%, also viel zu hoch für ein Random-Walk-Proposal.

Mit sehr hoher Varianz im Random Walk erhält man eine sehr niedrige Akzeptanzrate:

## Akzeptanzrate gesamt: 12

Tuning

• Daher tunt man den Algorithmus: Finde den optimalen Wert für \(c\)!
• Für Random Walk Proposal gelten Akzeptanzraten zwischen ca. 30% und 60% als optimal.
• Für andere Metropolis-Hastings-Schritte sind dagegen sehr hohe Akzeptanzraten oft in Ordnung, z.B. hier beim Ziehen von \(\tau\) durch eine Approximation der Full Conditional.

Hier tunen wir wie folgt:

• Algorithmus läuft für einige Iterationen (hier 50)

• Berechne die Akzeptanzrate

• Ist die Akzeptanzrate zu hoch: erhöhe Random-Walk-Varianz

• Ist die Akzeptanzrate zu niedrig: verringere Random-Walk-Varianz

• In beiden Fällen: starte Algorithmus erneut (übernimm aber aktuelle Werte der Parameter als Startwerte)

• Tuning und Burn-In überlagern sich dabei. Solange die Kette noch nicht “im richtigen Bereich” ist, ist ein anderer Tuning-Parameter eventuell effizienter.

• Auch andere Proposals können getunt werden.

Poisson-Lognormal mit Tuning

mcmc.tuning<-poisson.lognormal.mcmc(sumx,n,do.tuning=TRUE)

## Akzeptanzrate Tuning: 94

## Neuer Varianzparameter: 0.188

## Akzeptanzrate Tuning: 98

## Neuer Varianzparameter: 0.368

## Akzeptanzrate Tuning: 92

## Neuer Varianzparameter: 0.678

## Akzeptanzrate Tuning: 72

## Neuer Varianzparameter: 0.976

## Akzeptanzrate Tuning: 68

## Neuer Varianzparameter: 1.328

## Akzeptanzrate Tuning: 50

## Neuer Varianzparameter: 1.328

## Akzeptanzrate Tuning: 56

## Neuer Varianzparameter: 1.487

## Akzeptanzrate Tuning: 62

## Neuer Varianzparameter: 1.844

## Akzeptanzrate Tuning: 46

## Neuer Varianzparameter: 1.844

## Akzeptanzrate gesamt: 51.9

Diese Kette zeigt nun ein gutes Mischungsverhalten (good mixing).

Running mean plots

• Eine Möglichkeit, die Konvergenz der Kette bzw. des Posteriori-Erwartungswert zu visualisieren, ist der running mean.
• Wir plotten die Iteration gegen den Mittelwert des Parameters bis zu dieser Iteration
• der Mittelwert konvergiert nach dem Gesetz der großen Zahlen (bzw. dem Ergodensatz) gegen den Erwartungswert:

par(mfrow=c(1,2))
I=1000
plot(cumsum(mcmc.tuning[,1])/1:I,type="l", ylab="Mittelwert", xlab="Iteration")
plot(cumsum(mcmc.tuning[,2])/1:I,type="l", ylab="Mittelwert", xlab="Iteration")

Die Visualisierung ist dabei teilweise etwas irreführend. Eventuell handelt es sich nur um Bewegungen in den hinteren Nachkommastellen.

Auto correlation function

Den ACF-Plot kennne wir schon. coda hat dafür eine eigene Funktion. Hier der Vergleich des Algorithmus mit schlechtem und mit gutem Mixing:

acfplot(mcmc.simple)

acfplot(mcmc.tuning)

Gelman-Rubin-Diagnostik

• Idee: vergleiche die Varianz von mehreren parallel gelaufenen Ketten

• Bei Konvergenz sollte die Varianz in den Ketten der Varianz zwischen den Ketten entsprechen

• Alternativ ist auch die Aufteilung einer Kette in mehrere Teilketten möglich

• Within Chain Variance (unterschätzt die wahre Varianz, wenn Ketten noch nicht konvergiert):

\[ W=\frac{1}{m}\sum_{j=1}^ms_j^2, s_j^2=\frac{1}{n-1}\sum_{i=1}^n(\theta_{ij}-\bar{\theta_j})^2 \]

• Between Chain Variance

\[ B=\frac{n}{m-1}\sum_{j=1}^m(\bar{\theta}_j-\bar{\bar{\theta}})^2; \bar{\bar{\theta}}=\frac{1}{m}\sum_{j=1}^m\bar{\theta} = j \]

• Geschätzte Varianz

\[ \hat{Var}(\theta)=(1-\frac{1}{n}W+\frac{1}{n}B) \]

• Potential scale reduction factor

\[ \hat{R}=\sqrt\frac{\bar{Var}(\theta)}{W} \]

• Ist der Potential scale reduction factor \(\hat{R}\) zu groß (\(>1.1\)), sollten die Ketten länger laufen.
• Die Erweiterung ist der multivariate potential scale reduction factor, der alle Parameter gleichzeitig erfasst.

Anwendung

Zuerst als Beispiel unser MCMC-Algorithmus ohne Tuning. Wir lassen den Algorithmus fünf Mal laufen. Das Ergebnis wandeln wir in ein mcmc.list-Objekt um.

mc.list<-parallel::mclapply(rep(sumx, 5),
            poisson.lognormal.mcmc,n=n)
mc.list<-as.mcmc.list(mc.list)

Der Plot zeigt die fünf Ketten in fünf verschiedenen Farben. Die Dichte (rechts) wird dagegen aus den Ziehungen von allen fünf Ketten zusammen zusammen geschätzt.

• Der Traceplot zeigt zusätzlich eine durchgezogene geglätte Linie der Kette. Idealerweise sollte diese für alle Ketten konstant den gleichen Wert haben.
• Wir sehen bereits am Traceplot von \(\mu\), dass die Ketten hier teil unterschiedliche Ergebnisse liefern.
• Der Traceplot von \(\lambda\) sieht dagegen gut aus. VORSICHT! Den Ziehungen von \(\lambda\) ist trotzdem nicht zu trauen. Der Algorithmus musss insgesamt konvergieren, also für alle Parameter gut mischen!

gelman.diag

Die Funktion gelman.diag berechnet den potential scale reduction factor

gelman.diag(mc.list)

## Potential scale reduction factors:
## 
##        Point est. Upper C.I.
## mu           1.81       3.04
## lambda       1.02       1.03
## 
## Multivariate psrf
## 
## 1.79

• Der psrf (bzw. ) von \(\mu\) ist deutlich zu hoch.
• Auch der multivariate psrf ist zu hoch. Dieser ist insbesondere dann interessant, wenn man sehr viele Parameter hat.

Was macht man, wenn der psrf zu hoch ist?

• Mehr Iterationen (eventuell auch längerer Burn-In).
• Besserer Algorithmus.

Mehr Iterationen

mc.list<-parallel::mclapply(rep(sumx, 5),
            poisson.lognormal.mcmc,n=n,I=20000)
gelman.diag(mc.list)

## Potential scale reduction factors:
## 
##        Point est. Upper C.I.
## mu           1.02       1.06
## lambda       1.00       1.00
## 
## Multivariate psrf
## 
## 1.03

Algorithmus mit Tuning

Nur 1000 Iterationen (plus Iterationen für Tuning und burn-in):

mc.list<-parallel::mclapply(rep(sumx, 5),
            poisson.lognormal.mcmc,n=n, do.tuning=TRUE)
gelman.diag(mc.list)

## Potential scale reduction factors:
## 
##        Point est. Upper C.I.
## mu           1.01       1.03
## lambda       1.01       1.02
## 
## Multivariate psrf
## 
## 1.02

Geweke-Diagnostik

• Idee: Teste, ob zwei Teile einer Kette aus der selben Verteilung stammen (Teste Differenz der Mittelwerte)
• In der Regel erste 10% und letzte 50% der Kette
• geweke.diag: berechnet eine Z-Statistik für jeden Parameter

geweke.diag(mcmc.simple)

## 
## Fraction in 1st window = 0.1
## Fraction in 2nd window = 0.5 
## 
##      mu  lambda 
## -1.6617 -0.1108

pnorm(geweke.diag(mcmc.simple)$z)

##         mu     lambda 
## 0.04828307 0.45586997

Hier sind die ersten 10% und die letzten 50% ziemlich wahrscheinlich nicht aus der selben Verteilung.

geweke.diag(mcmc.tuning)

## 
## Fraction in 1st window = 0.1
## Fraction in 2nd window = 0.5 
## 
##       mu   lambda 
##  0.81908 -0.01463

pnorm(geweke.diag(mcmc.tuning)$z)

##        mu    lambda 
## 0.7936285 0.4941638

Hier vermutlich schon.

Heidelberg und Welch Diagnostik

Die Heidelberg-Welch-Diagnostik testet, ob die Kette lang genug ist, um den Mittelwert adäquat zu schätzen. Test wird für jeden Parameter separat ausgeführt.

• Teste, ob die Kette aus einer stationären Verteilung kommt
• Zuerst: Cramer-von Mises-Test mit Niveau \(\alpha\) auf ganzer Kette
• Falls Nullhypothese abgelehnt, verwirf erste 10%, teste erneut; falls abgelehnt, verwirf erste 20%, …, bis zu 50%
• Falls immer noch abgelehnt, muss die Kette länger sein.
• Falls Nullhypothese angenommen: Half-width-test.
• Berechne \((1-\alpha)\)-Kredibilitätsintervall (KI) für den Mittelwert (nur aus den letzen Iterationen, falls erste verworfen wurden).
• Teststatistik: Hälfte der Breite (Half-width) des Kredibilitätsintervall geteilt durch geschätzen Mittelwert

print(heidel.diag(mcmc.simple))

##                                      
##        Stationarity start     p-value
##        test         iteration        
## mu     passed       1         0.409  
## lambda passed       1         0.908  
##                                 
##        Halfwidth Mean  Halfwidth
##        test                     
## mu     failed     1.74 0.329    
## lambda passed    10.56 0.107

print(heidel.diag(mcmc.tuning))

##                                      
##        Stationarity start     p-value
##        test         iteration        
## mu     passed       1         0.271  
## lambda passed       1         0.993  
##                                
##        Halfwidth Mean Halfwidth
##        test                    
## mu     passed     2.3 0.1433   
## lambda passed    10.5 0.0962

• Beim Algorithmus ohne Tuning reicht die Anzahl der Iterationen nicht, die Teststatistik bei \(\mu\) ist zu groß
• Beim Algorithmus mit Tuning reicht die Anzahl der Iterationen zur Schätzung des Mittelwerts.

Raftery und Lewis Diagnostik

• Wir wollen ein Posteriori-Quantil \(q\) schätzen.
• Frage: Wieviele Iterationen \(N\) und welchen burn-in \(M\) brauchen wir, um mit Wahrscheinlichkeit \(s\) innerhalb der Toleranz \(r\) das Quantil \(q\) zu schätzen?
• Nach einer Pilotkette lassen wir die längere Kette entsprechend laufen.

raftery.diag schätzt die Länge des burn-in und die Anzahl der Iterationen aus der Pilotkette:

raftery.diag(mcmc.simple, 
                   q = 0.05, r = 0.02, s = 0.9)

## 
## Quantile (q) = 0.05
## Accuracy (r) = +/- 0.02
## Probability (s) = 0.9 
##                                               
##         Burn-in  Total Lower bound  Dependence
##         (M)      (N)   (Nmin)       factor (I)
##  mu     45       3612  322          11.20     
##  lambda 3        392   322           1.22

raftery.diag(mcmc.tuning, 
                   q = 0.05, r = 0.02, s = 0.9)

## 
## Quantile (q) = 0.05
## Accuracy (r) = +/- 0.02
## Probability (s) = 0.9 
##                                               
##         Burn-in  Total Lower bound  Dependence
##         (M)      (N)   (Nmin)       factor (I)
##  mu     15       1391  322          4.32      
##  lambda 3        361   322          1.12

• Für den getunten Algorithmus reicht die Anzahl der Iterationen schon fast für die gewünschte Genauigkeit.
• Für bessere Genauigkeit muss schon die Pilotkette länger sein.