apTreeshape (version 1.5-0.1)

mcmc_eta: Inference of the alpha and eta parameters

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.

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

# NOT RUN {
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)
# }
# NOT RUN {

# }