rmbayesgl {hwep}

Bayes test for random mating using genotype log-likelihoods

Description

Bayes test for random mating using genotype log-likelihoods

Usage

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().

Author(s)

David Gerard

References

• Gerard D (2022). "Bayesian tests for random mating in autopolyploids." bioRxiv. doi:10.1101/2022.08.11.503635.

Examples

## Not run: 
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")


## End(Not run)