Bayes test for random mating using genotype log-likelihoods

rmbayesgl(
  gl,
  method = c("stan", "gibbs"),
  lg = TRUE,
  alpha = NULL,
  beta = NULL,
  type = c("auto", "allo"),
  chains = 2,
  cores = 1,
  iter = 2000,
  warmup = floor(iter/2),
  ...
)

Arguments

gl: A matrix of genotype log-likelihoods. The rows index the individuals and the columns index the genotypes. So gl[i,k] is the genotype log-likelihood for individual i at dosage k-1. We assume the natural log is used (base e).
method: Should we use Stan ("stan") or Gibbs sampling ("gibbs")?
lg: A logical. Should we return the log Bayes factor (TRUE) or the Bayes factor (FALSE)?
alpha: The concentration hyperparameters of the gamete frequencies under the null of random mating. Should be length ploidy/2 + 1.
beta: The concentration hyperparameters of the genotype frequencies under the alternative of no random mating. Should be length ploidy + 1.
type: If alpha is NULL, then the default priors depend on if you have autopolyploids ("auto") or allopolyploids ("allo").
chains: Number of MCMC chains. Almost always 1 is enough, but we use 2 as a default to be conservative.
cores: Number of cores to use.
iter: Total number of iterations.
warmup: Number of those iterations used in the burnin.
...: Control arguments passed to sampling().

References

Gerard D (2023). "Bayesian tests for random mating in polyploids." Molecular Ecology Resources, In press. doi:10.1111/1755-0998.13856 .

Author

David Gerard

Examples


if (FALSE) {
set.seed(1)
ploidy <- 4

## Simulate under the null ----
p <- stats::runif(ploidy / 2 + 1)
p <- p / sum(p)
q <- stats::convolve(p, rev(p), type = "open")

## See BF increases
nvec <- c(stats::rmultinom(n = 1, size = 10, prob = q))
gl <- simgl(nvec = nvec)
rmbayesgl(gl = gl)
rmbayesgl(gl = gl, method = "gibbs")

nvec <- c(stats::rmultinom(n = 1, size = 100, prob = q))
gl <- simgl(nvec = nvec)
rmbayesgl(gl = gl)
rmbayesgl(gl = gl, method = "gibbs")

nvec <- c(stats::rmultinom(n = 1, size = 1000, prob = q))
gl <- simgl(nvec = nvec)
rmbayesgl(gl = gl)
rmbayesgl(gl = gl, method = "gibbs")

## Simulate under the alternative ----
q <- stats::runif(ploidy + 1)
q <- q / sum(q)

## See BF decreases
nvec <- c(stats::rmultinom(n = 1, size = 10, prob = q))
gl <- simgl(nvec = nvec)
rmbayesgl(gl = gl)
rmbayesgl(gl = gl, method = "gibbs")

nvec <- c(stats::rmultinom(n = 1, size = 100, prob = q))
gl <- simgl(nvec = nvec)
rmbayesgl(gl = gl)
rmbayesgl(gl = gl, method = "gibbs")

nvec <- c(stats::rmultinom(n = 1, size = 1000, prob = q))
gl <- simgl(nvec = nvec)
rmbayesgl(gl = gl)
rmbayesgl(gl = gl, method = "gibbs")

}