Description Usage Arguments Value Author(s) References Examples
Run the Bayesian inference of the clade age-richness index alpha and the clade abundance richness index eta
1 2 3 |
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) |
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.
Odile Maliet, Fanny Gascuel & Amaury Lambert
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
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 | 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)
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.