ewenss.calc: Direct calculation of Ewens's sampling formula with optional...

In wesleycrouse/TIMBR: Tree-based Inference of Multiallelism via Bayesian Regression

Description

Directly calculates probabilties for allelic series (partitions) under Ewen's sampling formula, optionally informed by a user-specified tree. Trees must be in coalescent units for appropriate inference.

Usage

ewenss.calc(tree, prior.alpha, stop.on.error = F)

Arguments

tree	either a user-specified tree of class "phylo" ("multiPhylo"), detailed in the 'ape' package, or an integer with the number of leaves to be partitioned
prior.alpha	prior type c("fixed","gamma","beta.prime") for the concentration parameter, see examples for format
stop.on.error	stop function if error is encountered when using 'integrate'. errors related to roundoff and small values may occur during edge cases. Note that function returns warning if error in total probability is >1 percent regardless of this setting

Value

list of allelic series IDs and probabilities, formatted as prior.M object for TIMBR function

Examples

#specifying hyperparameters for gamma prior using calc.concentration.prior
hyperparam <- calc.concentration.prior(8, 0.05, 0.01)
prior.alpha <- list(type="gamma", shape=hyperparam[1], rate=hyperparam[2])

#running the sampler without user-specified trees; compare with target prior probabilities
prior.M <- ewenss.calc(8, prior.alpha)
exp(prior.M$ln.probs[prior.M$M.IDs=="0,0,0,0,0,0,0,0"])
exp(prior.M$ln.probs[prior.M$M.IDs=="0,1,2,3,4,5,6,7"])

#running the sampler with a user-specified tree and fixed concentration parameter; compare with tree structure
tree <- ape::rcoal(8, LETTERS[1:8])
ape::plot.phylo(tree)
prior.alpha <- list(type="fixed", alpha=1)
prior.M <- ewenss.calc(tree, prior.alpha)
head(prior.M$M.IDs)
head(exp(prior.M$ln.probs))

wesleycrouse/TIMBR documentation built on Feb. 19, 2021, 7:31 a.m.