rmbayes: Bayes test for random mating with known genotypes

In hwep: Hardy-Weinberg Equilibrium in Polyploids

Bayes test for random mating with known genotypes

Description

Bayes test for random mating with known genotypes

Usage

rmbayes(
  nvec,
  lg = TRUE,
  alpha = NULL,
  beta = NULL,
  nburn = 10000,
  niter = 10000,
  type = c("auto", "allo")
)

Arguments

nvec	A vector containing the observed genotype counts, where nvec[[i]] is the number of individuals with genotype i-1. This should be of length ploidy+1.
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.
nburn	The number of iterations in the Gibbs sampler to burn-in.
niter	The number of sampling iterations in the Gibbs sampler.
type	If alpha is NULL, then the default priors depend on if you have autopolyploids ("auto") or allopolyploids ("allo").

Author(s)

David Gerard

References

• Gerard D (2022). "Bayesian tests for random mating in autopolyploids." bioRxiv. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1101/2022.08.11.503635")}.

Examples

set.seed(1)
ploidy <- 8

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

## See BF increase
nvec <- c(stats::rmultinom(n = 1, size = 100, prob = q))
rmbayes(nvec = nvec)

nvec <- c(stats::rmultinom(n = 1, size = 1000, prob = q))
rmbayes(nvec = nvec)

nvec <- c(stats::rmultinom(n = 1, size = 10000, prob = q))
rmbayes(nvec = nvec)

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

## See BF decrease
nvec <- c(stats::rmultinom(n = 1, size = 100, prob = q))
rmbayes(nvec = nvec)

nvec <- c(stats::rmultinom(n = 1, size = 1000, prob = q))
rmbayes(nvec = nvec)

nvec <- c(stats::rmultinom(n = 1, size = 10000, prob = q))
rmbayes(nvec = nvec)

hwep documentation built on May 31, 2023, 9:06 p.m.