gsi_mcmc_fb: MCMC from the fully Bayesian GSI model for pi and the...

In rubias: Bayesian Inference from the Conditional Genetic Stock Identification Model

Description

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.

Usage

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.

Value

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)

rubias documentation built on Feb. 10, 2022, 1:06 a.m.