fit_ClaDS: Fit ClaDS to a phylogeny

In hmorlon/PANDA: Phylogenetic ANalyses of DiversificAtion

Fit ClaDS to a phylogeny

Description

Performs the inferrence of branch-specific speciation rates and the model's hyper parameters for the model with constant extinction rate (ClaDS1) or constant turnover rate (ClaDS2).

Usage

fit_ClaDS(tree,sample_fraction,iterations, thin = 50, file_name = NULL, it_save = 1000,
                     model_id="ClaDS2", nCPU = 1, mcmcSampler = NULL, ...)

Arguments

tree	An object of class 'phylo'
sample_fraction	The sampling fraction for the clade on which the inference is performed.
iterations	Number of steps in the MCMC, should be a multiple of thin.
thin	Number of iterations between two chain state's recordings.
file_name	Name of the file in which the result will be saved. Use file_name = NULL (the default) to disable this option.
it_save	Number of iteration between each backup of the result in file_name.
model_id	"ClaDS1" for constant extinction rate, "ClaDS2" (the default) for constant turnover rate.
nCPU	The number of CPUs to use. Should be either 1 or 3.
mcmcSampler	Optional output of fit_ClaDS to continue an already started run.
...	Optional arguments, see details.

Details

This function uses a blocked Differential Evolution (DE) MCMC sampler, with sampling from the past of the chains (Ter Braak, 2006; ter Braak and Vrugt, 2008). This sampler is self-adaptive because proposals are generated from the past of the chains. In this sampler, three chains are run simultaneously. Block updates is implemented by first drawing the number of parameters to be updated from a truncated geometric distribution with mean 3, drawing uniformly which parameter to update, and then following the normal DE algorithm.

The available optional arguments are :

Nchain

Number of MCMC chains (default to 3).

res_ClaDS0

The output of ClaDS0 to use as a startpoint. If NULL (the default) a random startpoint is used for the branch-specific speciation rates for each chain.

l0

The starting value for lambda_0 (not used if res_ClaDS0 != NULL).

s0

The starting value for sigma (not used if res_ClaDS0 != NULL).

nlambda

Number of subdivisions for the rate space discretization (use in the likelihood computation). Default to 1000.

nt

Number of subdivisions for the time space discretization (use in the likelihood computation). Default to 30.

Value

A 'list' object with fields :

post	The posterior function.
startvalue	The starting value for the MCMC.
numPars	The number of parameter in the model, including the branch-specific speciation rates.
Nchain	The number of MCMC chains ran simultaneously.
currentLPs	The current values of the logposterior for th Nchains chains.
proposalGenerator	The proposal distribution for the MCMC sampler.
former	The last output of post for each of the chains.
thin	Number of iterations between two chain state's recordings.
alpha_effect	A vector of size nrow(tree$edge), where the ith element is the number of branches on the path from the crown of the tree and branch i (used internally in other functions).
consoleupdates	The frequency at which the sampler state should be printed.
likelihood	The likelihood function, used internally.
relToAbs	A function mapping the relative changes in speciation rates to the absolute speciation rates for the object phylo, used internally.

Author(s)

O. Maliet

References

Ter Braak, C. J. 2006. A markov chain monte carlo version of the genetic algorithm differential evolution: easy bayesian computing for real parameter spaces. Statistics and Computing 16:239- 249.

ter Braak, C. J. and J. A. Vrugt. 2008. Differential evolution markov chain with snooker updater and fewer chains. Statistics and Computing 18:435-446.

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

fit_ClaDS0, plot_ClaDS_chains.

Examples

if(test){
data("Caprimulgidae")

sample_fraction = 0.61

sampler = fit_ClaDS(Caprimulgidae, sample_fraction, 1000, thin = 50, 
          file_name = NULL, model_id="ClaDS2", nCPU = 1)
plot_ClaDS_chains(sampler)

# continue the same run 
sampler = fit_ClaDS(Caprimulgidae, sample_fraction, 50, mcmcSampler = sampler)




# plot the result of the analysis (saved in "Caprimulgidae_ClaDS2", after thinning)

data("Caprimulgidae_ClaDS2")

# plot the mcmc chains
plot_ClaDS_chains(Caprimulgidae_ClaDS2$sampler)

# extract the Maxima A Posteriori for each parameter
maps = getMAPS_ClaDS(Caprimulgidae_ClaDS2$sampler, thin = 1)
print(paste0("sigma = ", maps[1], " ; alpha = ", 
  maps[2], " ; epsilon = ", maps[3], " ; l_0 = ", maps[4] ))
  
# plot the infered branch specific speciation rates
plot_ClaDS_phylo(Caprimulgidae_ClaDS2$tree, maps[-(1:4)])
}

hmorlon/PANDA documentation built on April 24, 2024, 3:27 a.m.