fit_ClaDS0 function

Infer ClaDS0's parameter on a phylogeny

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

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.

Returns

A mcmc.list object with the three MCMC chains.

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

Author(s)

O. Maliet

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)])
     
}