dd_SR_ML {DDD} | R Documentation |
Maximization of the loglikelihood under a diversity-dependent diversification model with a shift in the parameters
Description
This function computes the maximum likelihood estimates of the parameters of a diversity-dependent diversification model with shifting parameters at time t = tshift for a given set of phylogenetic branching times. It also outputs the corresponding loglikelihood that can be used in model comparisons.
Usage
dd_SR_ML(
brts,
initparsopt = c(0.5, 0.1, 2 * (1 + length(brts) + missnumspec), 2 * (1 + length(brts) +
missnumspec), max(brts)/2),
parsfix = NULL,
idparsopt = c(1:3, 6:7),
idparsfix = NULL,
idparsnoshift = (1:7)[c(-idparsopt, (-1)^(length(idparsfix) != 0) * idparsfix)],
res = 10 * (1 + length(brts) + missnumspec),
ddmodel = 1,
missnumspec = 0,
cond = 1,
btorph = 1,
soc = 2,
allbp = FALSE,
tol = c(0.001, 1e-04, 1e-06),
maxiter = 1000 * round((1.25)^length(idparsopt)),
changeloglikifnoconv = FALSE,
optimmethod = "subplex",
num_cycles = 1,
methode = "analytical",
verbose = FALSE
)
Arguments
brts |
A set of branching times of a phylogeny, all positive |
initparsopt |
The initial values of the parameters that must be optimized |
parsfix |
The values of the parameters that should not be optimized |
idparsopt |
The ids of the parameters that must be optimized, e.g. 1:7
for all parameters. The ids are defined as follows: |
idparsfix |
The ids of the parameters that should not be optimized, e.g. c(1,3,4,6) if lambda and K should not be optimized, but only mu. In that case idparsopt must be c(2,5,7). The default is to fix all parameters not specified in idparsopt. |
idparsnoshift |
The ids of the parameters that should not shift; This can only apply to ids 4, 5 and 6, e.g. idparsnoshift = c(4,5) means that lambda and mu have the same values before and after tshift |
res |
sets the maximum number of species for which a probability must be computed, must be larger than 1 + length(brts) |
ddmodel |
sets the model of diversity-dependence: |
missnumspec |
The number of species that are in the clade but missing in the phylogeny |
cond |
Conditioning: |
btorph |
Sets whether the likelihood is for the branching times (0) or the phylogeny (1) |
soc |
Sets whether stem or crown age should be used (1 or 2) |
allbp |
Sets whether a search should be done with various initial conditions, with tshift at each of the branching points (TRUE/FALSE) |
tol |
Sets the tolerances in the optimization. Consists of: |
maxiter |
Sets the maximum number of iterations in the optimization |
changeloglikifnoconv |
if TRUE the loglik will be set to -Inf if ML does not converge |
optimmethod |
Method used in optimization of the likelihood. Current default is 'subplex'. Alternative is 'simplex' (default of previous versions) |
num_cycles |
the number of cycles of opimization. If set at Inf, it will do as many cycles as needed to meet the tolerance set for the target function. |
methode |
The method used to solve the master equation, default is 'analytical' which uses matrix exponentiation; alternatively numerical ODE solvers can be used, such as 'odeint::runge_kutta_cash_karp54'. These were used in the package before version 3.1. |
verbose |
Show the parameters and loglikelihood for every call to the loglik function |
Details
The output is a dataframe containing estimated parameters and maximum loglikelihood. The computed loglikelihood contains the factor q! m!/(q + m)! where q is the number of species in the phylogeny and m is the number of missing species, as explained in the supplementary material to Etienne et al. 2012.
Value
lambda_1 |
gives the maximum likelihood estimate of lambda before the shift |
mu_1 |
gives the maximum likelihood estimate of mu before the shift |
K_1 |
gives the maximum likelihood estimate of K before the shift |
lambda_2 |
gives the maximum likelihood estimate of lambda after the shift |
mu_2 |
gives the maximum likelihood estimate of mu after the shift |
K_2 |
gives the maximum likelihood estimate of K after the shift |
t_shift |
gives the time of the shift |
loglik |
gives the maximum loglikelihood |
df |
gives the number of estimated parameters, i.e. degrees of feedom |
conv |
gives a message on convergence of optimization; conv = 0 means convergence |
Note
The optimization may get trapped in local optima. Try different starting values to search for the global optimum.
Author(s)
Rampal S. Etienne & Bart Haegeman
References
- Etienne, R.S. et al. 2012, Proc. Roy. Soc. B 279: 1300-1309,
doi: 10.1098/rspb.2011.1439
- Etienne, R.S. & B. Haegeman 2012. Am. Nat.
180: E75-E89, doi: 10.1086/667574
See Also
dd_SR_loglik
, dd_ML
,
dd_KI_ML
,
Examples
cat("This will estimate parameters for a sets of branching times brts without conditioning.\n")
cat("The tolerance of the optimization is set ridiculously high to make runtime fast.\n")
cat("In real applications, use the default or more stringent settings for tol.\n")
brts = 1:10
dd_SR_ML(brts = brts, initparsopt = c(0.4581, 1E-6, 17.69, 11.09, 8.9999), idparsopt = c(1:3,6,7),
idparsfix = NULL, parsfix = NULL, idparsnoshift = c(4,5), cond = 0,
tol = c(1E-1,1E-1,1E-1),optimmethod = 'simplex'
)