fitMLE_bdK: Maximum likelihood fit of the birth-death k model on a...

In sophia-lambert/UDivEvo: A package for Unknown Diversity Diversification Analysis in Evolution

Maximum likelihood fit of the birth-death k model on a phylogeny

Description

Fits the birth-death k model (a constant-time homogeneous birth-death model under k-sampling) to a rooted ultrametric phylogeny using Maximum likelihood estimation. Precisely, it uses the Nelder-Mead algorithm to optimise the likelihood. The birth-death process is conditioned on the starting time of the process tot_time and the survival of the process at present time as well as having k extant sampled tips at present. This function can fit the birth-death k model on a stem or crown phylogeny. This function is specifically adapted for diversification analysis on phylogenies on which the sampling probability is unknown.

Usage

fitMLE_bdK(
  phylo,
  tot_time,
  r,
  epsi,
  cond = "crown",
  YULE = FALSE,
  dt = 0,
  rel.tol = 1e-10,
  tuned_dichotomy = TRUE,
  brk = 2000
)

Arguments

phylo	Object of class phylo. A rooted ultrametric phylogeny of class phylo. The rooted ultrametric phylogeny can have polytomie(s) (i.e. non binary tree).
tot_time	Numeric. The stem or crown age (also called MRCA) of the phylogeny depending on the conditioning of the process specified (see cond argument) and the phylogeny used accordingly. The stem age of the phylogeny can be computed using max(TreeSim::getx(phylo))+phylo$root.edge (note that the phylo$root.edge needs to be known) and the crown age of the phylogeny can be computed using max(TreeSim::getx(phylo)).
r	Numeric vector. The net diversification rate r starting value(s).
epsi	Numeric vector. The turnover rate ε starting value(s).
cond	Character. Specifying the conditioning of the birth-death process. Two conditioning are available, either cond = "crown" (the default option) if the phylogeny used is a crown phylogeny or cond = "stem" if the phylogeny used is a stem phylogeny.
YULE	Logical. If TRUE, the extinction rate μ thus the turnover rate ε are fixed to 0 and the net diversification rate r equals the speciation rate λ. If FALSE (the default option), the turnover rate ε is not fixed to 0 and is thus inferred.
dt	Numeric. If dt = 0, the integral on the sampling probability is computed using the R stats::integrate function. If dt≥0, the integral of the sampling probability is performed manually using a piece-wise constant approximation. dt represents the length of the interval on which the function integrated is assumed to be constant. For manual integral, advised value of dt are 1e-3 to 1e-5.
rel.tol	Numeric. This argument is only used if dt = 0. This represents the relative accuracy requested when the integral is performed using the stats::integrate function. Typically .Machine$double.eps^0.25 is used but a value of 1e-10 (the default value) has been tested and performs well.
tuned_dichotomy	Logical. If TRUE, when the log likelihood of the model is equal to non finite value due to approximations, a dichotomy search is performed to find a tuning parameter that will be used for getting a finite value of the log likelihood. If TRUE, the log likelihood will take longer to calculate. Else if FALSE, no dichotomy search is performed; if the log likelihood is equal to non finite value due to approximations, the log likelihood will take this non finite value for the corresponding parameters.
brk	Numeric. This argument is only used if tuned_dichotomy = TRUE. The number of steps used in the dichotomy search. Typically the value 200 is sufficient to avoid non finite values. In some case if the log likelihood is still equal to non finite value, the brk value 2000 will be required for more tuning but it will rarely take a larger value.

Details

This function will fit the birth-death k-sampling model and the function will infer the net diversification rate r and the turnover rate ε. Note that the starting values should be adapted to the number of parameters (here 2, r and ε). This function is specifically intended to be used on phylogenies with unknown or highly uncertain global diversity estimates (the sampling probability is not known with accuracy). Note that the sampling probability is never estimated and that this function is not able to evaluate negative rates.

Value

Returns an object of class MLE_bdK. This MLE_bdK object is a list containing the name of the birth-death model performed, the log likelihood of the data knowing the parameters, the akaike information criterion corrected and the inferred parameters.

Author(s)

Sophia Lambert

See Also

likelihood_bdK and fitMCMC_bdK

Examples

# Creating a phylogeny with 0.05 net diversification rate and 0.5 turnover rate.

set.seed(1234)
tree1 <- TESS::tess.sim.age(1, 100, lambda = 0.1, mu = 0.05, MRCA = TRUE, samplingProbability = 0.5)[[1]]
plot(tree1, root.edge = TRUE)

# Creating variables to give to arguments

tot_time <- max(TreeSim::getx(tree1))
Ntips <- ape::Ntip(tree1)
lamb_moments <- log(Ntips)/tot_time

# Fitting the birth-death k-sampling model

res_fitMCMC_M2 <- fitMLE_bdK(phylo = tree1,
                             tot_time = tot_time,
                             r = lamb_moments,
                             epsi = 0,
                             cond = "crown", YULE = FALSE,
                             dt = 0, rel.tol = 1e-10,
                             tuned_dichotomy = TRUE,
                             brk = 2000)

res_fitMLE_M2$parameters

sophia-lambert/UDivEvo documentation built on Sept. 27, 2022, 11:05 p.m.