fit_ClaDS0 {RPANDA}

Infer ClaDS0's parameter on a phylogeny

Description

Infer branch-specific speciation rates and the model's hyper parameters for the pure-birth model

Usage

fit_ClaDS0(tree, name, pamhLocalName = "pamhLocal",  
            iteration = 1e+07, thin = 20000, update = 1000, 
            adaptation = 10, seed = NULL, nCPU = 3)

Arguments

tree

An object of class 'phylo'.

name

The name of the file in which the results will be saved. Use name = NULL to disable this option.

pamhLocalName

The function is writing in a text file to make the execution quicker, this is the name of this file.

iteration

Number of iterations after which the gelman factor is computed and printed. The function stops if it is below 1.05

thin

Number of iterations between two chain state's recordings.

update

Number of iterations between two adjustments of the proposal parameters during the adaptation phase of the sampler.

adaptation

Number of times the proposal is adjusted during the adaptation phase of the sampler.

seed

An optional seed for the MCMC run.

nCPU

The number of CPUs to use. Should be either 1 or 3.

Details

This function uses a Metropolis within Gibbs MCMC sampler with a bactrian proposal (ref) with an initial adaptation phase. During this phase, the proposal is adjusted "adaptation" times every "update" iterations to reach a goal acceptance rate of 0.3.

To monitor convergence, 3 independant MCMC chains are run simultaneously and the Gelman statistics is computed every "iteration" iterations. The inference is stopped when the maximum of the one dimentional Gelman statistics (computed for each of the parameters) is below 1.05.

Value

A mcmc.list object with the three MCMC chains.

Author(s)

O. Maliet

References

Maliet O., Hartig F. and Morlon H. 2019, A model with many small shifts for estimating species-specific diversificaton rates, Nature Ecology and Evolution, doi 10.1038/s41559-019-0908-0

See Also

getMAPS_ClaDS0, plot_ClaDS0_chains, fit_ClaDS

Examples

set.seed(1)


if(test){

obj= sim_ClaDS( lambda_0=0.1,    
                mu_0=0.5,      
                sigma_lamb=0.7,         
                alpha_lamb=0.90,     
                condition="taxa",    
                taxa_stop = 20,    
                prune_extinct = TRUE)  

tree = obj$tree
speciation_rates = obj$lamb[obj$rates]
extinction_rates = obj$mu[obj$rates]

plot_ClaDS_phylo(tree,speciation_rates)

sampler = fit_ClaDS0(tree=tree,        
              name="ClaDS0_example.Rdata",      
              nCPU=1,               
              pamhLocalName = "local",
              iteration=500000,
              thin=2000,
              update=1000, adaptation=5) 
              
# extract the Maximum A Posteriori for each of the parameters
MAPS = getMAPS_ClaDS0(tree, sampler, thin = 10)

# plot the simulated (on the left) and inferred speciation rates (on the right)
# on the same color scale
plot_ClaDS_phylo(tree, speciation_rates, MAPS[-(1:3)])
     
}