mcmc_alpha: Inference of the alpha parameter

In apTreeshape: Analyses of Phylogenetic Treeshape

Description

Run the Bayesian inference of the clade age-richness index alpha

Usage

mcmc_alpha(tree, epsilon, beta, niter, ini = 0, V = 0.1, 
            chain = NULL, verbose = 10, silent = TRUE, Nadapt = Inf, 
            NadaptMin = 10, NadaptMax=Inf, ma = -4, Ma = 4, proposal = "bactrian", 
            accOpt = 0.3, Vmin = 0.001)

Arguments

tree	A phylo object
epsilon	Minimum size of unsampled splits (see appendix 1)
beta	Imbalance index
niter	Number of iterations in the mcmc
ini	Initial alpha value (default to 0)
V	Initial scaling value for the mcmc proposal (default to 0.1)
chain	Former mcmc chain (if NULL (the default), a new one is started)
verbose	Number of iterations after which the state of the mcmc is printed (if silent == FALSE)
silent	If TRUE (the default) the state of the mcmc is not printed
Nadapt	Number of iterations between each proposal scalling (default to Inf)
NadaptMin	Minimum nmber of iterations before the first proposal scalling (default to 10)
NadaptMax	Number of iterations after which the proposal stops being scalled (default to Inf)
ma	Minimal alpha value (default to -4)
Ma	Maximal alpha value (default to 4)
proposal	Shape of the proposal. Can be "bactrian" (the default, ref), "uniform", or "normal"
accOpt	Optimal acceptance value (default to 0.3)
Vmin	Minimal scaling value for the mcmc proposal (default to 0.001)

Value

A list with a mcmc field contening the resulting chain. The other fields are only used to resume runing the inference if the chain has to be completed.

Author(s)

Odile Maliet, Fanny Gascuel & Amaury Lambert

References

Maliet O., Gascuel F., Lambert A. (2018) Ranked tree shapes, non-random extinctions and the loss of phylogenetic diversity, bioRxiv 224295, doi: https://doi.org/10.1101/224295

Examples

ntip=30
set.seed(123)
tree=simulate_tree(epsilon = 0.01,alpha = -1,beta = 0,N = ntip,equal.ab = TRUE)
beta=maxlik.betasplit(tree,up=10)$max_lik
plot(tree)
  
niter=1000

## Not run: 
chain=mcmc_alpha(tree,epsilon=0.01,beta=beta,niter=600,V = c(0.1),ini=c(0),
                 verbose = 100,silent = FALSE,Nadapt = 100,NadaptMin = 100)
  
# Continue the same chain
chain=mcmc_alpha(tree,epsilon=0.01,beta=beta,niter=400,verbose = 100,silent = FALSE,
                 chain = chain,Nadapt = 100,NadaptMin = 500,NadaptMax = 700)
  
thinned=mcmc(chain$mcmc[seq(200,1000,10),])
plot(thinned)
da=density(thinned[,"alpha"])
MPa=da$x[which.max(da$y)]
print(MPa)


## End(Not run)

apTreeshape documentation built on Jan. 8, 2021, 2:07 a.m.