MCMC_BB  R Documentation 
This function initiates the Markov chain Monte Carlo (MCMC) for the betabinomial BEDASSLE model. The betabinomial model allows populations to diverge from the model's expectations based on their location and their neighbors.
MCMC_BB(counts, sample_sizes, D, E, k, loci, delta, aD_stp, aE_stp, a2_stp, phi_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 (D_{i,j}). This can be twodimensional Euclidean distance, or greatcircle distance, or, in fact, any positive definite matrix (deriving, for instance, from a resistance distance). However, note that the algorithm silently restricts the prior on the alpha parameters, and specifically the alpha_2 parameter, to the part of parameter space that results in valid covariance matrices; in the case of twodimensional Euclidean distances, this will not happen, since any value of alpha_2 between 0 and 2 is valid (see Guillot et al.'s "Valid covariance models for the analysis of geographical genetic variation" for more detail on this). 
E 
Pairwise ecological distance(s) (E_{i,j}), which may be continuous (e.g. 
difference in elevation) or binary (same or opposite side of some hypothesized
barrier to gene flow). Users may specify one or more ecological distance matrices.
If more than one is specified, they should be formatted as a 
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). 
phi_stp 
The scale of the tuning parameter on the phi parameters. 
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 betabinomial likelihood function may generate "NaNs produced" warnings, #so temporarily disable warnings. op < options("warn") options(warn = 1) #Call the Markov chain Monte Carlo for the overdispersion model ## Not run: MCMC_BB( 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, phi_stp = mcmc.operators$phi_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)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.