,
Bart Haegeman [aut] ,
Hanno Hildenbrandt [ctb] ,
Giovanni Laudanno [ctb] | Title: | Diversity-Dependent Diversification |
|---|---|
| Description: | Implements maximum likelihood and bootstrap methods based on the diversity-dependent birth-death process to test whether speciation or extinction are diversity-dependent, under various models including various types of key innovations. See Etienne et al. 2012, Proc. Roy. Soc. B 279: 1300-1309, <DOI:10.1098/rspb.2011.1439>, Etienne & Haegeman 2012, Am. Nat. 180: E75-E89, <DOI:10.1086/667574>, Etienne et al. 2016. Meth. Ecol. Evol. 7: 1092-1099, <DOI:10.1111/2041-210X.12565> and Laudanno et al. 2021. Syst. Biol. 70: 389–407, <DOI:10.1093/sysbio/syaa048>. Also contains functions to simulate the diversity-dependent process. |
| Authors: | Rampal S. Etienne [aut, cre] ,
Bart Haegeman [aut] ,
Hanno Hildenbrandt [ctb] ,
Giovanni Laudanno [ctb] |
| Maintainer: | Rampal S. Etienne <[email protected]> |
| License: | GPL-3 |
| Version: | 5.2.2 |
| Built: | 2024-10-27 04:19:31 UTC |
| Source: | https://github.com/rsetienne/ddd |
This function computes loglikelihood of a diversity-independent diversification model for a given set of branching times and parameter values.
bd_loglik( pars1, pars2, brts, missnumspec, methode = "odeint::runge_kutta_cash_karp54" )bd_loglik( pars1, pars2, brts, missnumspec, methode = "odeint::runge_kutta_cash_karp54" )
pars1 |
Vector of parameters:
|
pars2 |
Vector of model settings:
|
brts |
A set of branching times of a phylogeny, all positive |
missnumspec |
The number of species that are in the clade but missing in the phylogeny |
methode |
The method used to solve the master equation, default is 'odeint::runge_kutta_cash_karp54'. |
The loglikelihood
Rampal S. Etienne, Bart Haegeman & Cesar Martinez
- 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
bd_loglik(pars1 = c(0.5,0.1), pars2 = c(0,1,1,0,2), brts = 1:10, missnumspec = 0)bd_loglik(pars1 = c(0.5,0.1), pars2 = c(0,1,1,0,2), brts = 1:10, missnumspec = 0)
This function computes the maximum likelihood estimates of the parameters of a diversity-independent diversification model for a given set of phylogenetic branching times. It also outputs the corresponding loglikelihood that can be used in model comparisons.
bd_ML( brts, initparsopt = c(0.1, 0.05 * (tdmodel <= 1) + 10 * (length(brts) + missnumspec) * (tdmodel > 1)), idparsopt = c(1, 2 + (tdmodel > 1)), idparsfix = (1:4)[-idparsopt], parsfix = rep(0, 4)[idparsfix], missnumspec = 0, tdmodel = 0, cond = 1, btorph = 1, soc = 2, tol = c(0.001, 1e-04, 1e-06), maxiter = 1000 * round((1.25)^length(idparsopt)), changeloglikifnoconv = FALSE, optimmethod = "subplex", num_cycles = 1, methode = "odeint::runge_kutta_cash_karp54", verbose = FALSE )bd_ML( brts, initparsopt = c(0.1, 0.05 * (tdmodel <= 1) + 10 * (length(brts) + missnumspec) * (tdmodel > 1)), idparsopt = c(1, 2 + (tdmodel > 1)), idparsfix = (1:4)[-idparsopt], parsfix = rep(0, 4)[idparsfix], missnumspec = 0, tdmodel = 0, cond = 1, btorph = 1, soc = 2, tol = c(0.001, 1e-04, 1e-06), maxiter = 1000 * round((1.25)^length(idparsopt)), changeloglikifnoconv = FALSE, optimmethod = "subplex", num_cycles = 1, methode = "odeint::runge_kutta_cash_karp54", verbose = FALSE )
brts |
A set of branching times of a phylogeny, all positive |
initparsopt |
The initial values of the parameters that must be optimized |
idparsopt |
The ids of the parameters that must be optimized, e.g. 1:3
for intrinsic speciation rate, extinction rate and carrying capacity. The
ids are defined as follows: |
idparsfix |
The ids of the parameters that should not be optimized, e.g. c(1,3) if lambda0 and lambda1 should not be optimized, but only mu0 and mu1. In that case idparsopt must be c(2,4). The default is to fix all parameters not specified in idparsopt. |
parsfix |
The values of the parameters that should not be optimized |
missnumspec |
The number of species that are in the clade but missing in the phylogeny |
tdmodel |
Sets the model of time-dependence: |
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) |
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 under tdmodel = 4, default is 'odeint::runge_kutta_cash_karp54'. |
verbose |
Show the parameters and loglikelihood for every call to the loglik function |
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.
A dataframe with the following elements:
lambda0 |
gives the maximum likelihood estimate of lambda0 |
mu0 |
gives the maximum likelihood estimate of mu0 |
lambda1 |
gives the maximum likelihood estimate of lambda1 |
mu1 |
gives the maximum likelihood estimate of mu1 |
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 |
Rampal S. Etienne & Bart Haegeman
- 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
cat("Estimating parameters for a set of branching times brts with the default settings:") brts = 1:20 bd_ML(brts = brts, cond = 1)cat("Estimating parameters for a set of branching times brts with the default settings:") brts = 1:20 bd_ML(brts = brts, cond = 1)
Converting a set of branching times to a phylogeny
brts2phylo(times, root = FALSE, tip.label = NULL)brts2phylo(times, root = FALSE, tip.label = NULL)
times |
Set of branching times |
root |
When root is FALSE, the largest branching time will be assumed to be the crown age. When root is TRUE, it will be the stem age. |
tip.label |
Tip labels. If set to NULL, the labels will be t1, t2, etc. |
phy |
A phylogeny of the phylo type |
Rampal S. Etienne
- 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
Convolution of two vectors
conv(x, y)conv(x, y)
x |
first vector |
y |
second vector |
vector that is the convolution of x and y
Rampal S. Etienne
- 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
conv(1:10,1:10)conv(1:10,1:10)
This function computes loglikelihood of a diversity-dependent diversification model for a given set of branching times and parameter values where the diversity-dependent dynamics of a subclade decouple from the dynamics of the main clade at time t_d, potentially accompanied by a shift in parameters.
dd_KI_loglik( pars1, pars2, brtsM, brtsS, missnumspec, methode = "odeint::runge_kutta_cash_karp54" )dd_KI_loglik( pars1, pars2, brtsM, brtsS, missnumspec, methode = "odeint::runge_kutta_cash_karp54" )
pars1 |
Vector of parameters: |
pars2 |
Vector of model settings: |
brtsM |
A set of branching times of the main clade in the phylogeny, all positive |
brtsS |
A set of branching times of the subclade in the phylogeny, all positive |
missnumspec |
The number of species that are in the clade but missing in the phylogeny. One can specify the sum of the missing species in main clade and subclade or a vector c(missnumspec_M,missnumspec_S) with missing species in main clade and subclade respectively. |
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. |
The loglikelihood
Rampal S. Etienne & Bart Haegeman
- 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
dd_KI_ML, dd_loglik
dd_SR_loglik
pars1 = c(0.25,0.12,25.51,1.0,0.16,8.61,9.8) pars2 = c(200,1,0,18.8,1,2) missnumspec = 0 brtsM = c(25.2,24.6,24.0,22.5,21.7,20.4,19.9,19.7,18.8,17.1,15.8,11.8,9.7,8.9,5.7,5.2) brtsS = c(9.6,8.6,7.4,4.9,2.5) dd_KI_loglik(pars1,pars2,brtsM,brtsS,missnumspec)pars1 = c(0.25,0.12,25.51,1.0,0.16,8.61,9.8) pars2 = c(200,1,0,18.8,1,2) missnumspec = 0 brtsM = c(25.2,24.6,24.0,22.5,21.7,20.4,19.9,19.7,18.8,17.1,15.8,11.8,9.7,8.9,5.7,5.2) brtsS = c(9.6,8.6,7.4,4.9,2.5) dd_KI_loglik(pars1,pars2,brtsM,brtsS,missnumspec)
This function computes the maximum likelihood estimates of the parameters of a diversity-dependent diversification model with decoupling of the diversification dynamics of a subclade from the dynamics of the main clade for a given set of phylogenetic branching times of main clade and subclade and the time of splitting of the lineage that will form the subclade. It also outputs the corresponding loglikelihood that can be used in model comparisons.
dd_KI_ML( brtsM, brtsS, tsplit, initparsopt = c(0.5, 0.1, 2 * (1 + length(brtsM) + missnumspec[1]), 2 * (1 + length(brtsS) + missnumspec[length(missnumspec)]), (tsplit + max(brtsS))/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(c(brtsM, brtsS)) + sum(missnumspec)), ddmodel = 1, missnumspec = 0, cond = 1, soc = 2, tol = c(0.001, 1e-04, 1e-06), maxiter = 1000 * round((1.25)^length(idparsopt)), changeloglikifnoconv = FALSE, optimmethod = "subplex", num_cycles = 1, methode = "analytical", correction = TRUE, verbose = FALSE )dd_KI_ML( brtsM, brtsS, tsplit, initparsopt = c(0.5, 0.1, 2 * (1 + length(brtsM) + missnumspec[1]), 2 * (1 + length(brtsS) + missnumspec[length(missnumspec)]), (tsplit + max(brtsS))/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(c(brtsM, brtsS)) + sum(missnumspec)), ddmodel = 1, missnumspec = 0, cond = 1, soc = 2, tol = c(0.001, 1e-04, 1e-06), maxiter = 1000 * round((1.25)^length(idparsopt)), changeloglikifnoconv = FALSE, optimmethod = "subplex", num_cycles = 1, methode = "analytical", correction = TRUE, verbose = FALSE )
brtsM |
A set of branching times of the main clade in a phylogeny, all positive |
brtsS |
A set of branching times of the subclade in a phylogeny, all positive |
tsplit |
The branching time at which the lineage forming the subclade branches off, 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 + max(length(brtsM),length(brtsS)) |
ddmodel |
sets the model of diversity-dependence: |
missnumspec |
The number of species that are in the clade but missing in the phylogeny. One can specify the sum of the missing species in main clade and subclade or a vector c(missnumspec_M,missnumspec_S) with missing species in main clade and subclade respectively. |
cond |
Conditioning: |
soc |
Sets whether stem or crown age should be used (1 or 2); stem age only works when cond = 0 |
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. |
correction |
Sets whether the correction should be applied (TRUE) or not (FALSE) |
verbose |
Show the parameters and loglikelihood for every call to the loglik function |
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.
lambda_M |
gives the maximum likelihood estimate of lambda of the main clade |
mu_M |
gives the maximum likelihood estimate of mu of the main clade |
K_M |
gives the maximum likelihood estimate of K of the main clade |
lambda_2 |
gives the maximum likelihood estimate of lambda of the subclade |
mu_S |
gives the maximum likelihood estimate of mu of the subclade |
K_S |
gives the maximum likelihood estimate of K of the subclade |
t_d |
gives the time of the decoupling event |
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 |
The optimization may get trapped in local optima. Try different starting values to search for the global optimum.
Rampal S. Etienne & Bart Haegeman
- 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
dd_KI_loglik, dd_ML,
dd_SR_ML,
cat("This will estimate parameters for two sets of branching times brtsM, brtsS\n") cat("without conditioning.\n") cat("The tolerance of the optimization is set high so runtime is fast in this example.\n") cat("In real applications, use the default or more stringent settins for tol.\n") brtsM = 4:10 brtsS = seq(0.1,3.5,0.7) tsplit = 5 dd_KI_ML(brtsM = brtsM, brtsS = brtsS, tsplit = tsplit, idparsopt = c(1:3,6,7), initparsopt = c(0.885, 2e-14, 6.999, 6.848, 4.001), idparsfix = NULL, parsfix = NULL,idparsnoshift = c(4,5), cond = 0, tol = c(3E-1,3E-1,3E-1), optimmethod = 'simplex')cat("This will estimate parameters for two sets of branching times brtsM, brtsS\n") cat("without conditioning.\n") cat("The tolerance of the optimization is set high so runtime is fast in this example.\n") cat("In real applications, use the default or more stringent settins for tol.\n") brtsM = 4:10 brtsS = seq(0.1,3.5,0.7) tsplit = 5 dd_KI_ML(brtsM = brtsM, brtsS = brtsS, tsplit = tsplit, idparsopt = c(1:3,6,7), initparsopt = c(0.885, 2e-14, 6.999, 6.848, 4.001), idparsfix = NULL, parsfix = NULL,idparsnoshift = c(4,5), cond = 0, tol = c(3E-1,3E-1,3E-1), optimmethod = 'simplex')
Simulating a diversity-dependent diversification process where at a given time a new clade emerges with different inherent speciation rate and extinction rate and clade-level carrying capacity and with decoupled dynamics
dd_KI_sim(pars, age, ddmodel = 1)dd_KI_sim(pars, age, ddmodel = 1)
pars |
Vector of parameters: |
age |
Sets the crown age for the simulation |
ddmodel |
Sets the model of diversity-dependence: |
out |
A list with the following elements: The first element
is the tree of extant species in phylo format |
Rampal S. Etienne
- 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
dd_KI_sim(c(0.2,0.1,20,0.1,0.05,30,4),10)dd_KI_sim(c(0.2,0.1,20,0.1,0.05,30,4),10)
This function computes loglikelihood of a diversity-dependent diversification model for a given set of branching times and parameter values.
dd_loglik(pars1, pars2, brts, missnumspec, methode = "analytical")dd_loglik(pars1, pars2, brts, missnumspec, methode = "analytical")
pars1 |
Vector of parameters:
|
pars2 |
Vector of model settings:
|
brts |
A set of branching times of a phylogeny, all positive |
missnumspec |
The number of species that are in the clade but missing in the phylogeny |
methode |
The method used to solve the master equation, default is 'analytical' which uses matrix exponentiation; alternatively numerical ODE solvers can be used. Before version 3.1 these were solvers from the deSolve package such as 'lsoda' and 'ode45'. Currently solvers from odeint are used, such as 'odeint::runge_kutta_cash_karp54', 'odeint::runge_kutta_fehlberg78', 'odeint::runge_kutta_dopri5', or odeint::bulirsch_stoer'. The first two are recommended in most cases. |
The loglikelihood
Rampal S. Etienne & Bart Haegeman
- 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
dd_ML, dd_SR_loglik,
dd_KI_loglik
dd_loglik(pars1 = c(0.5,0.1,100), pars2 = c(100,1,1,1,0,2), brts = 1:10, missnumspec = 0)dd_loglik(pars1 = c(0.5,0.1,100), pars2 = c(100,1,1,1,0,2), brts = 1:10, missnumspec = 0)
This function computes the maximum likelihood and the associated estimates of the parameters of a diversity-dependent diversification model for a given set of phylogenetic branching times. It then performs a bootstrap likelihood ratio test of the diversity-dependent (DD) model against the constant-rates (CR) birth-death model. Finally, it computes the power of this test.
dd_LR( brts, initparsoptDD, initparsoptCR, missnumspec, outputfilename = NULL, seed = 42, endmc = 1000, alpha = 0.05, plotit = TRUE, res = 10 * (1 + length(brts) + missnumspec), ddmodel = 1, cond = 1, btorph = 1, soc = 2, tol = c(0.001, 1e-04, 1e-06), maxiter = 2000, changeloglikifnoconv = FALSE, optimmethod = "subplex", methode = "analytical" )dd_LR( brts, initparsoptDD, initparsoptCR, missnumspec, outputfilename = NULL, seed = 42, endmc = 1000, alpha = 0.05, plotit = TRUE, res = 10 * (1 + length(brts) + missnumspec), ddmodel = 1, cond = 1, btorph = 1, soc = 2, tol = c(0.001, 1e-04, 1e-06), maxiter = 2000, changeloglikifnoconv = FALSE, optimmethod = "subplex", methode = "analytical" )
brts |
A set of branching times of a phylogeny, all positive |
initparsoptDD |
The initial values of the parameters that must be optimized for the diversity-dependent (DD) model: lambda_0, mu and K |
initparsoptCR |
The initial values of the parameters that must be optimized for the constant-rates (CR) model: lambda and mu |
missnumspec |
The number of species that are in the clade but missing in the phylogeny |
outputfilename |
The name (and location) of the file where the output will be saved. Default is no save. |
seed |
The seed for the pseudo random number generator for simulating the bootstrap data |
endmc |
The number of bootstraps |
alpha |
The significance level of the test |
plotit |
Boolean to plot results or not |
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: |
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) |
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) |
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. |
The output is a list with 3 elements:
treeCR |
a list of trees generated under the constant-rates model using the ML parameters under the CR model |
treeDD |
a list of trees generated under the diversity-dependent model using the ML parameters under the diversity-dependent model |
out |
a dataframe with the
parameter estimates and maximum likelihoods for diversity-dependent and
constant-rates models |
pvalue |
p-value of the test |
LRalpha |
Likelihood ratio at the signifiance level alpha |
poweroftest |
power of the test for significance level alpha |
Rampal S. Etienne & Bart Haegeman
- Etienne, R.S. et al. 2016. Meth. Ecol. Evol. 7: 1092-1099,
doi: 10.1111/2041-210X.12565
- 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
This function computes the maximum likelihood estimates of the parameters of a diversity-dependent diversification model for a given set of phylogenetic branching times. It also outputs the corresponding loglikelihood that can be used in model comparisons.
dd_ML( brts, initparsopt = initparsoptdefault(ddmodel, brts, missnumspec), idparsopt = 1:length(initparsopt), idparsfix = (1:(3 + (ddmodel == 5)))[-idparsopt], parsfix = parsfixdefault(ddmodel, brts, missnumspec, idparsopt), res = 10 * (1 + length(brts) + missnumspec), ddmodel = 1, missnumspec = 0, cond = 1, btorph = 1, soc = 2, 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 )dd_ML( brts, initparsopt = initparsoptdefault(ddmodel, brts, missnumspec), idparsopt = 1:length(initparsopt), idparsfix = (1:(3 + (ddmodel == 5)))[-idparsopt], parsfix = parsfixdefault(ddmodel, brts, missnumspec, idparsopt), res = 10 * (1 + length(brts) + missnumspec), ddmodel = 1, missnumspec = 0, cond = 1, btorph = 1, soc = 2, 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 )
brts |
A set of branching times of a phylogeny, all positive |
initparsopt |
The initial values of the parameters that must be optimized |
idparsopt |
The ids of the parameters that must be optimized, e.g. 1:3
for intrinsic speciation rate, extinction rate and carrying capacity. The
ids are defined as follows: |
idparsfix |
The ids of the parameters that should not be optimized, e.g. c(1,3) if lambda and K should not be optimized, but only mu. In that case idparsopt must be 2. The default is to fix all parameters not specified in idparsopt. |
parsfix |
The values of the parameters that should not be optimized |
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) |
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 |
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.
lambda |
gives the maximum likelihood estimate of lambda |
mu |
gives the maximum likelihood estimate of mu |
K |
gives the maximum likelihood estimate of K |
r |
(only if ddmodel == 5) gives the ratio of linear dependencies in speciation and extinction rates |
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 |
Rampal S. Etienne & Bart Haegeman
- 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
dd_loglik, dd_SR_ML,
dd_KI_ML,
cat("Estimating the intrinsic speciation rate lambda and the carrying capacity K") cat("for a fixed extinction rate of 0.1, conditioning on clade survival and two missing species:") brts = 1:5 dd_ML(brts = brts,initparsopt = c(1.3078,7.4188), idparsopt = c(1,3), parsfix = 0.1, cond = 1, missnumspec = 2, tol = c(1E-3,1E-3,1E-4), optimmethod = 'simplex')cat("Estimating the intrinsic speciation rate lambda and the carrying capacity K") cat("for a fixed extinction rate of 0.1, conditioning on clade survival and two missing species:") brts = 1:5 dd_ML(brts = brts,initparsopt = c(1.3078,7.4188), idparsopt = c(1,3), parsfix = 0.1, cond = 1, missnumspec = 2, tol = c(1E-3,1E-3,1E-4), optimmethod = 'simplex')
This function computes the loglikelihood of a diversity-dependent diversification model for a given set of branching times and parameter values where the diversity-dependent dynamics of an innovative subclade have different parameters from the dynamics of the main clade from time t_d, but both are governed by the same carrying capacity and experience each other's diversity.
dd_MS_loglik( pars1, pars2, brtsM, brtsS, missnumspec, methode = "odeint::runge_kutta_cash_karp54" )dd_MS_loglik( pars1, pars2, brtsM, brtsS, missnumspec, methode = "odeint::runge_kutta_cash_karp54" )
pars1 |
Vector of parameters: |
pars2 |
Vector of model settings: |
brtsM |
A set of branching times of the main clade in the phylogeny, all positive |
brtsS |
A set of branching times of the subclade in the phylogeny, all positive |
missnumspec |
The number of species that are in the clade but missing in the phylogeny. One can specify the sum of the missing species in main clade and subclade or a vector c(missnumspec_M,missnumspec_S) with missing species in main clade and subclade respectively. |
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. |
The loglikelihood
Rampal S. Etienne & Bart Haegeman
- 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
dd_MS_ML, dd_loglik,
dd_KI_loglik, dd_SR_loglik
pars1 = c(0.2,0.1,40,1.0,0.1,9.8) pars2 = c(200,1,0,18.8,1,2) missnumspec = 0 brtsM = c(25.2,24.6,24.0,22.5,21.7,20.4,19.9,19.7,18.8,17.1,15.8,11.8,9.7,8.9,5.7,5.2) brtsS = c(9.6,8.6,7.4,4.9,2.5) dd_MS_loglik(pars1,pars2,brtsM,brtsS,missnumspec)pars1 = c(0.2,0.1,40,1.0,0.1,9.8) pars2 = c(200,1,0,18.8,1,2) missnumspec = 0 brtsM = c(25.2,24.6,24.0,22.5,21.7,20.4,19.9,19.7,18.8,17.1,15.8,11.8,9.7,8.9,5.7,5.2) brtsS = c(9.6,8.6,7.4,4.9,2.5) dd_MS_loglik(pars1,pars2,brtsM,brtsS,missnumspec)
This function computes the maximum likelihood estimates of the parameters of a diversity-dependent diversification model where the diversity-dependent dynamics of an innovative subclade have different parameters from the dynamics of the main clade from time t_d, but both are governed by the same carrying capacity and experience each other's diversity. Required isa given set of phylogenetic branching times of main clade and subclade and the time of splitting of the lineage that will form the subclade. The function also outputs the corresponding loglikelihood that can be used in model comparisons.
dd_MS_ML( brtsM, brtsS, tsplit, initparsopt = c(0.5, 0.1, 2 * (1 + length(brtsM) + length(brtsS) + sum(missnumspec)), (tsplit + max(brtsS))/2), parsfix = NULL, idparsopt = c(1:3, 6), idparsfix = NULL, idparsnoshift = (1:6)[c(-idparsopt, (-1)^(length(idparsfix) != 0) * idparsfix)], res = 10 * (1 + length(c(brtsM, brtsS)) + sum(missnumspec)), ddmodel = 1.3, missnumspec = 0, cond = 0, soc = 2, tol = c(0.001, 1e-04, 1e-06), maxiter = 1000 * round((1.25)^length(idparsopt)), changeloglikifnoconv = FALSE, optimmethod = "subplex", num_cycles = 1, methode = "ode45", correction = FALSE, verbose = FALSE )dd_MS_ML( brtsM, brtsS, tsplit, initparsopt = c(0.5, 0.1, 2 * (1 + length(brtsM) + length(brtsS) + sum(missnumspec)), (tsplit + max(brtsS))/2), parsfix = NULL, idparsopt = c(1:3, 6), idparsfix = NULL, idparsnoshift = (1:6)[c(-idparsopt, (-1)^(length(idparsfix) != 0) * idparsfix)], res = 10 * (1 + length(c(brtsM, brtsS)) + sum(missnumspec)), ddmodel = 1.3, missnumspec = 0, cond = 0, soc = 2, tol = c(0.001, 1e-04, 1e-06), maxiter = 1000 * round((1.25)^length(idparsopt)), changeloglikifnoconv = FALSE, optimmethod = "subplex", num_cycles = 1, methode = "ode45", correction = FALSE, verbose = FALSE )
brtsM |
A set of branching times of the main clade in a phylogeny, all positive |
brtsS |
A set of branching times of the subclade in a phylogeny, all positive |
tsplit |
The branching time at which the lineage forming the subclade branches off, 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 + max(length(brtsM),length(brtsS)) |
ddmodel |
sets the model of diversity-dependence: |
missnumspec |
The number of species that are in the clade but missing in the phylogeny. One can specify the sum of the missing species in main clade and subclade or a vector c(missnumspec_M,missnumspec_S) with missing species in main clade and subclade respectively. |
cond |
Conditioning: |
soc |
Sets whether stem or crown age should be used (1 or 2); stem age only works when cond = 0 |
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 in the ode solver. This can be either 'analytical' for explicit matrix exponentation or any of the solvers in the deSolve package. |
correction |
Sets whether the correction should be applied (TRUE) or not (FALSE) |
verbose |
Show the parameters and loglikelihood for every call to the loglik function |
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.
lambda_M |
gives the maximum likelihood estimate of lambda of the main clade |
mu_M |
gives the maximum likelihood estimate of mu of the main clade |
K_M |
gives the maximum likelihood estimate of K of the main clade |
lambda_2 |
gives the maximum likelihood estimate of lambda of the subclade |
mu_S |
gives the maximum likelihood estimate of mu of the subclade |
t_d |
gives the time of the key innovation event |
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 |
The optimization may get trapped in local optima. Try different starting values to search for the global optimum.
Rampal S. Etienne & Bart Haegeman
- 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
dd_MS_loglik, dd_ML,
dd_KI_ML, dd_SR_ML,
cat("This will estimate parameters for two sets of branching times brtsM, brtsS\n") cat("without conditioning.\n") cat("The tolerance of the optimization is set high so runtime is fast in this example.\n") cat("In real applications, use the default or more stringent settins for tol.\n") brtsM = 4:10 brtsS = seq(0.1,3.5,0.7) tsplit = 5 dd_MS_ML(brtsM = brtsM, brtsS = brtsS, tsplit = tsplit, idparsopt = c(1:3,6), initparsopt = c(0.885, 2e-14, 10, 4.001), idparsfix = NULL, parsfix = NULL, idparsnoshift = c(4,5), cond = 0, tol = c(3E-1,3E-1,3E-1))cat("This will estimate parameters for two sets of branching times brtsM, brtsS\n") cat("without conditioning.\n") cat("The tolerance of the optimization is set high so runtime is fast in this example.\n") cat("In real applications, use the default or more stringent settins for tol.\n") brtsM = 4:10 brtsS = seq(0.1,3.5,0.7) tsplit = 5 dd_MS_ML(brtsM = brtsM, brtsS = brtsS, tsplit = tsplit, idparsopt = c(1:3,6), initparsopt = c(0.885, 2e-14, 10, 4.001), idparsfix = NULL, parsfix = NULL, idparsnoshift = c(4,5), cond = 0, tol = c(3E-1,3E-1,3E-1))
Simulating a diversity-dependent diversification process where at a given time a new clade emerges with different inherent speciation rate and extinction rate
dd_MS_sim(pars, age, ddmodel = 1.3)dd_MS_sim(pars, age, ddmodel = 1.3)
pars |
Vector of parameters: |
age |
Sets the crown age for the simulation |
ddmodel |
Sets the model of diversity-dependence: |
out |
A list with the following elements: The first element
is the tree of extant species in phylo format |
Rampal S. Etienne
- 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
dd_MS_sim(c(0.2,0.1,20,0.1,0.05,4),10)dd_MS_sim(c(0.2,0.1,20,0.1,0.05,4),10)
This function computes loglikelihood of a diversity-dependent diversification model for a given set of branching times and parameter values where the diversity-dependent dynamics of subclades decouple from the dynamics of main clades, potentially accompanied by a shift in parameters.
dd_multiple_KI_loglik( pars1_list, pars2, brts_k_list, missnumspec_list, reltol = 1e-14, abstol = 1e-16, methode = "odeint::runge_kutta_cash_karp54" )dd_multiple_KI_loglik( pars1_list, pars2, brts_k_list, missnumspec_list, reltol = 1e-14, abstol = 1e-16, methode = "odeint::runge_kutta_cash_karp54" )
pars1_list |
list of paramater sets one for each rate regime (subclade). The parameters are: lambda (speciation rate), mu (extinction rate), and K (clade-level carrying capacity). |
pars2 |
Vector of model settings: |
brts_k_list |
list of matrices, one for each rate regime (subclade). Each matrix has in the first row the branching times including the shift/decoupling time and the present time (0) in negative time (i.e. 10 mya = -10). In the second row it has the number of lineages, i.e. starting at 2 for a phylogeny with a crown and increasing by one at each branching time and decreasing by one at each decoupling/shift time. The last element is the same as the second last. |
missnumspec_list |
list containing the number of missing species for each clade. If only a single number m of missing species is known for the entire phylogeny, then each element of the list should be 0:m. One can also create this from m using the function create_missnumspec_list |
reltol |
relative tolerance in integration of the ODE system, default at 1e-14 |
abstol |
tolerance tolerance in integration of the ODE system, default at 1e-16 |
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. |
Simulating the diversity-dependent diversification process
dd_sim(pars, age, ddmodel = 1)dd_sim(pars, age, ddmodel = 1)
pars |
Vector of parameters: |
age |
Sets the crown age for the simulation |
ddmodel |
Sets the model of diversity-dependence: |
out |
A list with the following four elements: The first
element is the tree of extant species in phylo format |
Rampal S. Etienne
- 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
dd_sim(c(0.2,0.1,20),10)dd_sim(c(0.2,0.1,20),10)
This function computes loglikelihood of a diversity-dependent diversification model for a given set of branching times and parameter values where the parameters are allowed to shift at time t = tshift
dd_SR_loglik(pars1, pars2, brts, missnumspec, methode = "analytical")dd_SR_loglik(pars1, pars2, brts, missnumspec, methode = "analytical")
pars1 |
Vector of parameters: |
pars2 |
Vector of model settings: |
brts |
A set of branching times of a phylogeny, all positive |
missnumspec |
The number of species that are in the clade but missing in the phylogeny |
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. |
The loglikelihood
Rampal S. Etienne & Bart Haegeman
- 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
dd_SR_ML, dd_loglik,
dd_KI_loglik
dd_SR_loglik(pars1 = c(0.2,0.1,50,0.2,0.1,70,5), pars2 = c(100,1,1,1,0,2), brts = 1:10, missnumspec = 0)dd_SR_loglik(pars1 = c(0.2,0.1,50,0.2,0.1,70,5), pars2 = c(100,1,1,1,0,2), brts = 1:10, missnumspec = 0)
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.
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 )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 )
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 |
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.
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 |
The optimization may get trapped in local optima. Try different starting values to search for the global optimum.
Rampal S. Etienne & Bart Haegeman
- 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
dd_SR_loglik, dd_ML,
dd_KI_ML,
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' )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' )
Simulating the diversity-dependent diversification process with a parameter shift at a certain time
dd_SR_sim(pars, age, ddmodel = 1)dd_SR_sim(pars, age, ddmodel = 1)
pars |
Vector of parameters: |
age |
Sets the crown age for the simulation |
ddmodel |
Sets the model of diversity-dependence: |
out |
A list with the following four elements: The first
element is the tree of extant species in phylo format |
Rampal S. Etienne
- 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
dd_SR_sim(c(0.2,0.1,20,0.2,0.1,40,5),10)dd_SR_sim(c(0.2,0.1,20,0.2,0.1,40,5),10)
Converting a table with speciation and extinction events to a set of branching times
L2brts(L, dropextinct = T)L2brts(L, dropextinct = T)
L |
Matrix of events as produced by dd_sim: |
dropextinct |
Sets whether the phylogeny should drop species that are extinct at the present |
brts |
A set of branching times |
Rampal S. Etienne
- 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
sim = dd_sim(c(0.2,0.1,20),10) phy = L2brts(sim$L) plot(phy)sim = dd_sim(c(0.2,0.1,20),10) phy = L2brts(sim$L) plot(phy)
Converting a table with speciation and extinction events to a phylogeny
L2phylo(L, dropextinct = T)L2phylo(L, dropextinct = T)
L |
Matrix of events as produced by dd_sim: |
dropextinct |
Sets whether the phylogeny should drop species that are extinct at the present |
phy |
A phylogeny of the phylo type |
Rampal S. Etienne
- 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
sim = dd_sim(c(0.2,0.1,20),10) phy = L2phylo(sim$L) plot(phy)sim = dd_sim(c(0.2,0.1,20),10) phy = L2phylo(sim$L) plot(phy)
A wrapper to use several optimization routines, currently only 'simplex' (a method adopted from Matlab, or 'subplex', from the R package subplex). The function is called from several packages by the same author.
optimizer( optimmethod = "simplex", optimpars = c(1e-04, 1e-04, 1e-06, 1000), num_cycles = 1, fun, trparsopt, jitter = 0, ... )optimizer( optimmethod = "simplex", optimpars = c(1e-04, 1e-04, 1e-06, 1000), num_cycles = 1, fun, trparsopt, jitter = 0, ... )
optimmethod |
The method to use for optimization, either 'simplex' or 'subplex' |
optimpars |
Parameters of the optimization: relative tolerance in function arguments, relative tolerance in function value, absolute tolerance in function arguments, and maximum number of iterations |
num_cycles |
Number of cycles of the optimization. When set to Inf, the optimization will be repeated until the result is, within the tolerance, equal to the starting values, with a maximum of 10 cycles. |
fun |
Function to be optimized |
trparsopt |
Initial guess of the parameters to be optimized |
jitter |
Perturbation of an initial parameter value when precisely equal to 0.5; this is only relevant when subplex is chosen. The default value is 0, so no jitter is applied. A recommended value when using it is 1E-5. |
... |
Any other arguments of the function to be optimimzed, or settings of the optimization routine |
out |
A list containing optimal function arguments
( |
.
Rampal S. Etienne
cat("No examples")cat("No examples")
Converting a phylogeny to a table with speciation and extinction events
phylo2L(phy)phylo2L(phy)
phy |
A phylogeny of the phylo type |
L |
Matrix of events as produced by dd_sim: |
Liang Xu
- 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
sim = dd_sim(c(0.2,0.1,20),10) phy = sim$tas L = phylo2L(phy) phy2 = L2phylo(L, dropextinct = FALSE) graphics::par(mfrow = c(1,3)) graphics::plot(phy) graphics::plot(phy2) graphics::plot(L2phylo(sim$L, dropextinct = FALSE))sim = dd_sim(c(0.2,0.1,20),10) phy = sim$tas L = phylo2L(phy) phy2 = L2phylo(L, dropextinct = FALSE) graphics::par(mfrow = c(1,3)) graphics::plot(phy) graphics::plot(phy2) graphics::plot(L2phylo(sim$L, dropextinct = FALSE))
Sampling in which zero probabilities are removed
rng_respecting_sample(x, size, replace, prob)rng_respecting_sample(x, size, replace, prob)
x |
either a vector of one or more elements from which to choose, or a positive integer. See ‘Details.’ |
size |
a non-negative integer giving the number of items to choose. |
replace |
should sampling be with replacement? |
prob |
a vector of probability weights for obtaining the elements of the vector being sampled. |
a vector of length size with elements drawn from either x or from the integers 1:x.
thanks to Pedro Neves for finding this feature in base::sample
Richel J.C. Bilderbeek
See sample for more details
# Number of draws n <- 1000 # Do normal sampling set.seed(42) draws_1 <- DDD:::rng_respecting_sample( 1:3, size = n, replace = TRUE, prob = c(1.0, 1.0, 1.0) ) # Do a sampling with one element of probability zero set.seed(42) draws_2 <- DDD:::rng_respecting_sample( 1:4, size = n, replace = TRUE, prob = c(1.0, 1.0, 1.0, 0.0) ) testit::assert(sum(draws_2 == 4) == 0) testit::assert(draws_1 == draws_2) # Use base sampling will give different results, # as it results in different RNG values set.seed(42) draws_3 <- sample( 1:4, size = n, replace = TRUE, prob = c(1.0, 1.0, 1.0, 0.0) ) testit::assert(sum(draws_3 == 4) == 0) testit::assert(!all(draws_1 == draws_3))# Number of draws n <- 1000 # Do normal sampling set.seed(42) draws_1 <- DDD:::rng_respecting_sample( 1:3, size = n, replace = TRUE, prob = c(1.0, 1.0, 1.0) ) # Do a sampling with one element of probability zero set.seed(42) draws_2 <- DDD:::rng_respecting_sample( 1:4, size = n, replace = TRUE, prob = c(1.0, 1.0, 1.0, 0.0) ) testit::assert(sum(draws_2 == 4) == 0) testit::assert(draws_1 == draws_2) # Use base sampling will give different results, # as it results in different RNG values set.seed(42) draws_3 <- sample( 1:4, size = n, replace = TRUE, prob = c(1.0, 1.0, 1.0, 0.0) ) testit::assert(sum(draws_3 == 4) == 0) testit::assert(!all(draws_1 == draws_3))
The standard round function in R rounds x.5 to the nearest even integer. This is odd behavior that is corrected in roundn
roundn(x, digits = 0)roundn(x, digits = 0)
x |
Number to be rounded |
digits |
Sets the number of decimals in rounding. |
n |
A number |
Rampal S. Etienne
round(2.5) roundn(2.5) round(3.5) roundn(3.5) round(2.65,digits = 1) roundn(2.65,digits = 1) round(2.75,digits = 1) roundn(2.75,digits = 1)round(2.5) roundn(2.5) round(3.5) roundn(3.5) round(2.65,digits = 1) roundn(2.65,digits = 1) round(2.75,digits = 1) roundn(2.75,digits = 1)
The standard sample function in R samples from n numbers when x = n. This is unwanted behavior when the size of the vector to sample from changes dynamically. This is corrected in sample2
sample2(x, size, replace = FALSE, prob = NULL)sample2(x, size, replace = FALSE, prob = NULL)
x |
A vector of one or more elements |
size |
A non-negative integer giving the number of items to choose. |
replace |
Should sampling be with replacement? |
prob |
A vector of probability weights for obtaining the elements of the vector being sampled. |
sam |
A vector of length |
Rampal S. Etienne
sample(x = 10,size = 5,replace = TRUE) sample2(x = 10,size = 5,replace = TRUE)sample(x = 10,size = 5,replace = TRUE) sample2(x = 10,size = 5,replace = TRUE)
Function to optimize target function using a simplex method adopted from Matlab
simplex(fun, trparsopt, optimpars, ...)simplex(fun, trparsopt, optimpars, ...)
fun |
Function to be optimized |
trparsopt |
Initial guess of the parameters to be optimized |
optimpars |
Parameters of the optimization: relative tolerance in function arguments, relative tolerance in function value, absolute tolerance in function arguments, and maximum number of iterations |
... |
Any other arguments of the function to be optimimzed, or settings of the optimization routine |
out |
A list containing optimal function arguments
( |
.
Rampal S. Etienne
cat("No examples")cat("No examples")
Simulates a phylogenetic tree branching according to a time-dependent process calibrated on the expected number of species under a diversity-dependent process over time.
td_sim(pars, age, ddmodel = 1, methode = "ode45")td_sim(pars, age, ddmodel = 1, methode = "ode45")
pars |
Vector of parameters:
|
age |
crown age of the tree to simulate, i.e. the simulation time. |
ddmodel |
the diversity-dependent model used as reference for the time-dependent model. |
methode |
The method used to solve the master equation.
See |
A list with the following four elements: The first element is the
tree of extant species in phylo format
The second element is the tree of all species, including extinct species,
in phylo format
The third element is a matrix of all species where
- the first column is the time at which a species is born
- the second column is the label of the parent of the species;
positive and negative values only indicate whether the species belongs to
the left or right crown lineage
- the third column is the label of the daughter species itself;
positive and negative values only indicate whether the species belongs to
the left or right crown lineage
- the fourth column is the time of extinction of the species. If this
equals -1, then the species is still extant.
César Martinez, Rampal S. Etienne
Function to transform pars in a way that is more useful for optimization: trpars <- sign(pars) * pars/(sign(pars) + pars);
transform_pars(pars)transform_pars(pars)
pars |
Parameters to be transformed |
Transformed parameters
Rampal S. Etienne
Function to untransform pars after optimization: pars <- sign(trpars) * trpars/(sign(trpars) - trpars);
untransform_pars(trpars)untransform_pars(trpars)
trpars |
Parameters to be untransformed |
Untransformed parameters
Rampal S. Etienne