Here are all the objects saved in the output:

library(BEDASSLE)
show(load("sim/CGW_4_MCMC_output1.Robj"))
##  [1] "last.params"   "LnL_thetas"    "LnL_counts"    "LnL"          
##  [5] "Prob"          "a0"            "aD"            "aE"           
##  [9] "a2"            "beta"          "samplefreq"    "ngen"         
## [13] "a0_moves"      "aD_moves"      "aE_moves"      "a2_moves"     
## [17] "thetas_moves"  "mu_moves"      "beta_moves"    "aD_accept"    
## [21] "aE_accept"     "a2_accept"     "thetas_accept" "mu_accept"    
## [25] "aD_stp"        "aE_stp"        "a2_stp"        "thetas_stp"   
## [29] "mu_stp"
fig.dim <- 5
require(knitr)
## Loading required package: knitr
opts_chunk$set(fig.height=fig.dim,fig.width=2*fig.dim,fig.align='center')

Diagnostics

First, let’s look at traces of the \(\alpha\) parameters and their acceptance rates:

layout(t(1:2))
plot(aD, xlab="MCMC generation", ylab="value", main="aD")
plot(aD_accept/aD_moves, xlab="MCMC generation", ylab="", main="aD acceptance", ylim=c(0,1))

plot(as.vector(aE), xlab="MCMC generation", ylab="value", main="aE")
plot(aE_accept/aE_moves, xlab="MCMC generation", ylab="", main="aE acceptance", ylim=c(0,1))

plot(a2, xlab="MCMC generation", ylab="value", main="a2")
plot(a2_accept/a2_moves, xlab="MCMC generation", ylab="", main="a2 acceptance", ylim=c(0,1))

And, the log-likelihood trace:

plot(-Prob, xlab="MCMC generation", main="log likelihood")

and mean acceptance rates for the allele-specific global mean (\(\mu\)) and population mean (\(\Theta\)) parameters:

plot(mu_accept/mu_moves, xlab="MCMC generation", ylab="", main="mu acceptance", ylim=c(0,1) )

plot(thetas_accept/thetas_moves, xlab="MCMC generation", ylab="", main="thetas acceptance", ylim=c(0,1) )

Here are autocorrelation traces of the \(\alpha\) parameters:

acf(a0,lag.max=1000,xlab='',ylab='')

acf(a2,lag.max=1000,xlab='',ylab='')

acf(as.vector(aE),lag.max=1000,xlab='',ylab='')

acf(aD,lag.max=1000,xlab='',ylab='')

Also, joint traces of parameters:

layout(matrix(1:6,nrow=3))
cols <- adjustcolor(rainbow(64)[cut(seq_along(a0),breaks=64)],0.5)
varnames <- c("aD","aE","a2","a0")
for (i in 1:3) {
    for (j in (i+1):4) {
        plot( as.vector(get(varnames[i])), as.vector(get(varnames[j])), col=cols, pch=20, xlab=varnames[i], ylab=varnames[j] ) 
    }
}

Results

Here’s the main thing we are interested in: the posterior distribution of the ratio of the relative strengths of \(\alpha_E\) and \(\alpha_D\):

hist(aE/aD,breaks=100,main="posterior of aE/aD ratio")