Example BEDASSLE run: HGDP

Here are all the objects saved in the output:

library(BEDASSLE)
show(load("hgdp-ex/him_MCMC_output.Robj"))

##  [1] "a0"            "aD"            "aE"            "a2"           
##  [5] "beta"          "a0_moves"      "aD_moves"      "aE_moves"     
##  [9] "a2_moves"      "beta_moves"    "mu_moves"      "thetas_moves" 
## [13] "aD_accept"     "aE_accept"     "a2_accept"     "mu_accept"    
## [17] "thetas_accept" "LnL"           "LnL_thetas"    "LnL_counts"   
## [21] "Prob"          "ngen"

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[-(1:1000)], xlab="MCMC generation", ylab="value", main="aD")
plot((aD_accept/aD_moves)[-(1:1000)], xlab="MCMC generation", ylab="", main="aD acceptance", ylim=c(0,1))

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

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

And, the log-likelihood trace:

plot(Prob[-(1:1000)], 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)[-(1:1000)], xlab="MCMC generation", ylab="", main="mu acceptance", ylim=c(0,1) )

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

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

acf(a0[-(1:1000)],lag.max=1000,xlab='',ylab='')

acf(a2[-(1:1000)],lag.max=1000,xlab='',ylab='')

acf(as.vector(aE)[-(1:1000)],lag.max=1000,xlab='',ylab='')

acf(aD[-(1:1000)],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]))[-(1:1000)], as.vector(get(varnames[j]))[-(1:1000)], 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)[-(1:1000)],breaks=100,main="posterior of aE/aD ratio")