fit_bd_in_past | R Documentation |
Fits the birth-death model with potentially time-varying rates and potentially missing extant species to a phylogeny, by maximum likelihood while excluding the recent past. Notations follow Morlon et al. PNAS 2011.
fit_bd_in_past(phylo, tot_time, time_stop, f.lamb, f.mu, lamb_par, mu_par, desc, tot_desc,
meth = "Nelder-Mead", cst.lamb = FALSE, cst.mu = FALSE,
expo.lamb = FALSE, expo.mu = FALSE, fix.mu = FALSE,
dt=0, cond = "crown")
phylo |
an object of type 'phylo' (see ape documentation) that does not include any recent speciation (i.e. no speciation events between time_stop and the present). |
time_stop |
the age of the phylogeny where to stop the birth-death process: it excludes the recent past (between the present and time_stop), while conditioning on the survival of the lineages from time_stop to the present. |
tot_time |
the age of the phylogeny (crown age, or stem age if known). If working with crown ages, tot_time is given by max(node.age(phylo)$ages). |
f.lamb |
a function specifying the hypothesized functional form (e.g. constant, linear, exponential, etc.) of the variation of the speciation rate |
f.mu |
a function specifying the hypothesized functional form (e.g. constant, linear, exponential, etc.) of the variation of the extinction rate |
lamb_par |
a numeric vector of initial values for the parameters of f.lamb to be estimated (these values are used by the optimization algorithm). The length of this vector is used to compute the total number of parameters in the model, so to fit a model with constant speciation rate (for example), lamb_par should be a vector of length 1. Otherwise aic values will be wrong. |
mu_par |
a numeric vector of initial values for the parameters of f.mu to be estimated (these values are used by the optimization algorithm). The length of this vector is used to compute the total number of parameters in the model, so to fit a model without extinction (for example), mu_par should be empty (vector of length 0). Otherwise aic values will be wrong. |
desc |
the number of lineages present at present in the reconstructed phylogenetic tree. |
tot_desc |
the total number of extant species (including in the unsampled ones). |
meth |
optimization to use to maximize the likelihood function, see optim for more details. |
cst.lamb |
logical: should be set to TRUE only if f.lamb is constant (i.e. does not depend on time) to use analytical instead of numerical computation in order to reduce computation time. |
cst.mu |
logical: should be set to TRUE only if f.mu is constant (i.e. does not depend on time) to use analytical instead of numerical computation in order to reduce computation time. |
expo.lamb |
logical: should be set to TRUE only if f.lamb is exponential to use analytical instead of numerical computation in order to reduce computation time. |
expo.mu |
logical: should be set to TRUE only if f.mu is exponential to use analytical instead of numerical computation in order to reduce computation time. |
fix.mu |
logical: if set to TRUE, the extinction rate |
dt |
the default value is 0. In this case, integrals in the likelihood are computed using R "integrate" function, which can be quite slow. If a positive dt is given as argument, integrals are computed using a piece-wise contant approximation, and dt represents the length of the intervals on which functions are assumed to be constant. For an exponential dependency of the speciation rate with time, we found that dt=1e-3 gives a good trade-off between precision and computation time. |
cond |
conditioning to use to fit the model:
|
The lengths of lamb_par and mu_par are used to compute the total number of parameters in the model, so to fit a model with constant speciation rate (for example), lamb_par should be a vector of length 1. Otherwise aic values will be wrong. In the f.lamb and f.mu functions, the first argument (time) runs from the present to the past. Hence, if the parameter controlling the variation of \lambda
with time is estimated to be positive (for example), this means that the speciation rate decreases from past to present. Note that abs(f.lamb) and abs(f.mu) are used in the likelihood computation as speciation and extinction rates should always be positive. A consequence of
this is that negative speciation/extinction rates estimates can be returned. They should be interpreted in aboslute terms. See Morlon et al. 2020 for a more detailed explanation.
a list with the following components
model |
the name of the fitted model |
LH |
the maximum log-likelihood value |
aicc |
the second order Akaike's Information Criterion |
lamb_par |
a numeric vector of estimated f.lamb parameters, in the same order as defined in f.lamb |
mu_par |
a numeric vector of estimated f.mu parameters, in the same order as defined in f.mu (if fix.mu is FALSE) |
H Morlon, E Lewitus, B Perez-Lamarque
Morlon, H., Parsons, T.L. and Plotkin, J.B. (2011) Reconciling molecular phylogenies with the fossil record Proc Nat Acad Sci 108: 16327-16332
Lewitus, E., Bittner, L., Malviya, S., Bowler, C., & Morlon, H. (2018) Clade-specific diversification dynamics of marine diatoms since the Jurassic Nature Ecology and Evolution, 2(11), 1715–1723
Perez-Lamarque, B., Öpik, M., Maliet, O., Afonso Silva, A., Selosse, M-A., Martos, F., Morlon, H. (2022) Analysing diversification dynamics using barcoding data: The case of an obligate mycorrhizal symbiont, Molecular Ecology 31: 3496–3512
fit_env_in_past
,fit_bd
,plot_fit_bd
, plot_dtt
library(ape)
library(phytools)
data(Cetacea)
plot(Cetacea)
tot_time<-max(node.age(Cetacea)$ages)
# slice the Cetaceae tree 10 Myr ago:
time_stop=10
sliced_tree <- Cetacea
sliced_sub_trees <- treeSlice(sliced_tree,slice = tot_time-time_stop, trivial=TRUE)
for (i in 1:length(sliced_sub_trees)){
if (Ntip(sliced_sub_trees[[i]])>1){
sliced_tree <- drop.tip(sliced_tree,
tip=sliced_sub_trees[[i]]$tip.label[2:Ntip(sliced_sub_trees[[i]])])
}}
for (i in which(node.depth.edgelength(sliced_tree)>(tot_time-time_stop))){
temp = sliced_tree$edge.length[which(sliced_tree$edge[,2]==i)]-time_stop
sliced_tree$edge.length[which(sliced_tree$edge[,2]==i)] <- temp
}
Ntip(sliced_tree) # 27 lineages present 10 Myr have survived until today
# Now we can fit birth-death models excluding the 10 last Myr
# Fit the pure birth model (no extinction) with a constant speciation rate
f.lamb <-function(t,y){y[1]}
f.mu<-function(t,y){0}
lamb_par<-c(0.09)
mu_par<-c()
result_cst <- fit_bd_in_past(sliced_tree, tot_time, time_stop, f.lamb, f.mu,
desc=Ntip(Cetacea), tot_desc=89, lamb_par, mu_par,
cst.lamb = TRUE, fix.mu=TRUE, dt=1e-3)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.