gsi_mcmc_fb function

MCMC from the fully Bayesian GSI model for pi and the individual posterior probabilities

Given a list of key parameters from a genetic dataset, this function samples values of pi and the posteriors for all the individuals. Each MCMC iteration includes a recalculation of the scaled genotype likelihood matrix, with baseline allele frequencies updated based on the previous iteration's allocations. It returns the output in a list.

gsi_mcmc_fb(
  par_list,
  Pi_init,
  lambda,
  reps,
  burn_in,
  sample_int_Pi,
  sample_int_PofZ
)

Arguments

• par_list: genetic data converted to the param_list format by tcf2param_list
• Pi_init: Starting value for the pi (collection mixture proportion) vector.
• lambda: the prior to be added to the collection allocations, in order to generate pseudo-count Dirichlet parameters for the simulation of a new pi vector
• reps: total number of reps (sweeps) to do.
• burn_in: how many reps to discard in the beginning when doing the mean calculation. They will still be returned in the traces if desired
• sample_int_Pi: the number of reps between samples being taken for Pi traces. If 0 no trace samples are taken
• sample_int_PofZ: the number of reps between samples being taken for the traces of posterior of each individual's origin. If 0 no trace samples are taken.

Returns

gsi_mcmc_fb returns a list of three. $mean lists the posterior means for collection proportions pi, for the individual posterior probabilities of assignment PofZ, and for the allele frequencies theta. $sd returns the posterior standard deviations for the same values.

If the corresponding sample_int variables are not 0, $trace contains samples taken from the Markov chain at intervals of sample_int_(variable) steps.

Examples

# this demonstrates it with scaled likelihoods computed from
# assignment of the reference samples

# we have to get the ploidies to pass to tcf2param_list
locnames <- names(alewife)[-(1:16)][c(TRUE, FALSE)]
ploidies <- rep(2, length(locnames))
names(ploidies) <- locnames

params <- tcf2param_list(alewife, 17, ploidies = ploidies)
lambda <- rep(1/params$C, params$C)
# use very short run and burn in so it doesn't take too long
# when checking on CRAN
mcmc <- gsi_mcmc_fb(params, lambda, lambda, 20, 5, 4, 4)