Assessing MCMCs

The first step in any MCMC is to determine whether it has been run for long enough. One way that we do this is to look at the probability or likelihood of each state sampled during the MCMC run. If the MCMC reaches some value and then fails to increase for a very long time this is consistent with it having run long enough. Often times we will also run the MCMC multiple times and make sure that the parameters being sampled after the MCMC has reached this point are largely the same between independent MCMC runs.

Lets read in some MCMC results

dat1 <- read.csv("birdmcmc1.csv")
dat2 <- read.csv("birdmcmc2.csv")
dat3 <- read.csv("birdmcmc3.csv")

Now with our data read in lets go ahead and assess each MCMC

plot(dat1$p, type="l")

This one looks like it reached stationarity very quickly and the vast majority of the MCMC run was sampling from the posterior.

plot(dat2$p, type="l")

Same thing for this one very quickly reached stationarity and then continued to sample from the posterior.

plot(dat3$p, type="l")

This one got stuck at a lower are of parameter space for a while before begining to sample from the posterior around iteration 1500.

Estimating parameters

To learn what our estimate of a parameter is we have to throw away the burnin portion of the MCMC and base our conclusion on only the posterior. For instance if we wanted to know the mean of asc1 and the desc1 rate in the first dataset we would do this as shown here:

post <- dat1[5001:10000,]
mean(post$asc1)
## [1] 46.39713
mean(post$desc1)
## [1] 45.59186

The impact of correctly throwing out the burnin will vary from analysis to analysis lets look at dataset 3 for an example of where it has a large impact. In this dataset lets look at ascending 1 (fission) rates with and without properly discarding the burnin

# if we throw away the burnin we get
mean(dat3$asc1[2001:4000])
## [1] 0.1204623
# if we faile to throw away the burnin we get
mean(dat3$asc1)
## [1] 0.1952375

Comparing parameters

Often times with Bayesian analyses we will answer our underlying questions by structuring them as comparisons between estimated parameters. Let do this with dataset 1. Specifically lets ask whether asc1 and asc2 are different from one another (these are two rates of chromosome fission). Because these are Bayesian we need to use credible intervals rather than confidence intervals lets use the package coda to use the HPDInterval function.

HPDinterval(as.mcmc(dat1$asc1[5001:10000]))
##         lower    upper
## var1 36.38109 56.60195
## attr(,"Probability")
## [1] 0.95
HPDinterval(as.mcmc(dat1$asc2[5001:10000]))
##         lower    upper
## var1 1.358938 3.902688
## attr(,"Probability")
## [1] 0.95

This shows us that the upper limit of the credible interval for asc2 is far below the lower limit of the credible interval of asc1. This means that we would say these are significantly different from each other. Below is an example of the type of plot that I would typcially make for somehting like this.

asc1 <- HPDinterval(as.mcmc(dat1$asc1[5001:10000]))[1:2]
asc2 <- HPDinterval(as.mcmc(dat1$asc2[5001:10000]))[1:2]
plot(density(dat1$asc1[5001:10000]),main="", xlab="fission rate", xlim=c(0,70),
     ylim=c(-0.05,.7))
lines(density(dat1$asc2[5001:10000]))
polygon(density(dat1$asc1[5001:10000]),col=rgb(1,0,0,.5))
polygon(density(dat1$asc2[5001:10000]),col=rgb(0,1,0,.5))
lines(x=asc1,y=c(-.025,-.025),lwd=3,col=rgb(1,0,0,.5))
lines(x=asc2,y=c(-.025,-.025),lwd=3,col=rgb(0,1,0,.5))
points(x=c(50,50),y=c(.65,.6), col=c(rgb(1,0,0,.5),rgb(0,1,0,.5)),pch=15)
text(x=c(50,50),y=c(.65,.6), c("State 1", "State 2"), pos=4, cex=.7)