mcmc_eta: Inference of the alpha and eta parameters

In apTreeshape: Analyses of Phylogenetic Treeshape

Description

Run the Bayesian inference of the clade age-richness index alpha and the clade abundance richness index eta

Usage

mcmc_eta(tree, epsilon, beta, ini = c(0, 1), V = c(0.1, 0.1), chain = NULL, niter, 
          verbose = 10, silent = TRUE, Nadapt = Inf, NadaptMin = 10, NadaptMax=Inf, 
          ma = -4, Ma = 4, me = 0.1, Me = 10, proposal = "bactrian", accOpt = 0.3)

Arguments

tree	A phylo object
epsilon	Minimum size of unsampled splits (see appendix 1)
beta	Imbalance index
ini	Initial alpha and eta values (default to c(0,1))
V	Initial scaling value for the mcmc proposal (default to c(0.1,0.1))
chain	Former mcmc chain (if NULL (the default), a new one is started)
niter	Number of iterations in the mcmc
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)
me	Minimal eta value (default to 0.1)
Me	Maximal eta value (default to 10)
proposal	Shape of the proposal. Can be "bactrian" (the default, ref), "uniform", or "normal"
accOpt	Optimal acceptance value (default to 0.3)

Value

Returns 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

seed=123
set.seed(seed)
ntip=30
tree=simulate_tree(epsilon = 0.001,alpha = -1,beta = 0,N = ntip,equal.ab = FALSE,eta =1.5)
beta=maxlik.betasplit(tree,up=10)$max_lik
extinctions = rank(tree$tip.ab)
tree$tip.label = rep(".", length(tree$tip.label))
plot.phylo(tree, show.node.label=TRUE, 
            cex=order(extinctions, seq(1,(tree$Nnode+1)))/
            ((tree$Nnode+1)/6), adj=0.1)

## Not run: 
chain=mcmc_eta(tree,epsilon=0.001,beta=beta,V = c(0.1,0.1),niter=600,ini=c(0,1),
                 verbose = 100,silent = FALSE,Nadapt = 100,NadaptMin = 100)
# The initialisation of the mcmc is quiet long because 
# we begin by drawing many unsampled intervals. 
# When this is done it gets quicker. 


chain=mcmc_eta(tree,epsilon=0.001,beta=beta,niter=400,verbose = 200,silent = FALSE,
                 chain = chain,Nadapt = 100,NadaptMin = 100,NadaptMax = 700)

thinned=mcmc(chain$mcmc[seq(200,1000,10),])
plot(thinned)
da=density(thinned[,"alpha"])
MPa=da$x[which.max(da$y)]
de=density(log(thinned[,"eta"]))
MPe=exp(de$x[which.max(de$y)])
print(MPa)
print(MPe)
## End(Not run)

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