MCMC {BEDASSLE}  R Documentation 
This function initiates the Markov chain Monte Carlo (MCMC) for the binomial BEDASSLE model.
MCMC(counts, sample_sizes, D, E, k, loci, delta, aD_stp, aE_stp, a2_stp, thetas_stp,
mu_stp, ngen, printfreq, savefreq, samplefreq, directory = NULL, prefix = "",
continue = FALSE, continuing.params = NULL)
counts 
A matrix of allelic count data, for which 
sample_sizes 
A matrix of sample sizes, for which 
D 
Pairwise geographic distance ( 
E 
Pairwise ecological distance(s) ( 
k 
The number of populations in the analysis. This should be equal to

loci 
The number of loci in the analysis. This should be equal to

delta 
The size of the "delta shift" on the offdiagonal elements of the parametric
covariance matrix, used to ensure its positivedefiniteness (even, for example,
when there are separate populations sampled at the same geographic/ecological
coordinates). This value must be large enough that the covariance matrix is
positivedefinite, but, if possible, should be smaller than the smallest off
diagonal distance elements, lest it have an undue impact on inference. If the
user is concerned that the delta shift is too large relative to the pairwise
distance elements in 
aD_stp 
The scale of the tuning parameter on aD (alphaD). The scale of the tuning
parameter is the standard deviation of the normal distribution from which small
perturbations are made to those parameters updated via a randomwalk sampler.
A larger value of the scale of the tuning parameter will lead to, on average,
larger proposed moves and lower acceptance rates (for more on acceptance rates,
see 
aE_stp 
The scale of the tuning parameter on aE (alphaE). If there are multiple ecological distances included in the analysis, there will be multiple alphaE parameters (one for each matrix in the list of E). These may be updated all with the same scale of a tuning parameter, or they can each get their own, in which case aE_stp should be a vector of length equal to the number of ecological distance variables. 
a2_stp 
The scale of the tuning parameter on a2 (alpha_2). 
thetas_stp 
The scale of the tuning parameter on the theta parameters. 
mu_stp 
The scale of the tuning parameter on mu. 
ngen 
The number of generations over which to run the MCMC (one parameter is updated at random per generation, with mu, theta, and phi all counting, for the purposes of updates, as one parameter). 
printfreq 
The frequency with which MCMC progress is printed to the screen. If

savefreq 
The frequency with which the MCMC saves its output as an R object ( 
samplefreq 
The thinning of the MCMC chain ( 
directory 
If specified, this points to a directory into which output will be saved. 
prefix 
If specified, this prefix will be added to all output file names. 
continue 
If 
continuing.params 
The list of parameter values used to initiate the MCMC if 
This function saves an MCMC output object at intervals specified by savefreq
.
This object may be ported into R working memory using the base function
load
.
As with any MCMC method, it is very important here to perform MCMC diagnosis and
evaluate chain mixing. I have provided a number of MCMC diagnosis graphing functions
for user convenience in visually assessing MCMC output. These include
plot_all_trace
;plot_all_marginals
;
plot_all_joint_marginals
;
and plot_all_acceptance_rates
. To evaluate model adequacy, users should
use posterior.predictive.sample
and
plot_posterior_predictive_sample
. These MCMC diagnosis/model adequacy
functions all call the standard MCMC output R object that the BEDASSLE MCMC generates
as their principal argument.
If users wish to start another MCMC run from where the current run left off, they
should use make.continuing.params
, and initiate the new run with option
continue = TRUE
and the continuing.params
list from the previous run
specified.
Gideon Bradburd
Bradburd, G.S., Ralph, P.L., and Coop, G.M. Disentangling the effects of geographic and ecological isolation on genetic differentiation. Evolution 2013.
#With the HGDP dataset and mcmc operators
data(HGDP.bedassle.data)
data(mcmc.operators)
#The value of delta may set off warnings,
#so temporarily disable warnings.
op < options("warn")
options(warn = 1)
#Call the Markov chain Monte Carlo for the standard model
## Not run:
MCMC(
counts = HGDP.bedassle.data$allele.counts,
sample_sizes = HGDP.bedassle.data$sample.sizes,
D = HGDP.bedassle.data$GeoDistance,
E = HGDP.bedassle.data$EcoDistance,
k = HGDP.bedassle.data$number.of.populations,
loci = HGDP.bedassle.data$number.of.loci,
delta = mcmc.operators$delta,
aD_stp = mcmc.operators$aD_stp,
aE_stp = mcmc.operators$aE_stp,
a2_stp = mcmc.operators$a2_stp,
thetas_stp = mcmc.operators$thetas_stp,
mu_stp = mcmc.operators$mu_stp,
ngen = mcmc.operators$ngen,
printfreq = mcmc.operators$printfreq,
savefreq = mcmc.operators$savefreq,
samplefreq = mcmc.operators$samplefreq,
directory = NULL,
prefix = mcmc.operators$prefix,
continue = FALSE,
continuing.params = NULL
)
## End(Not run)
#Reenable warnings
options(op)