simulate_diversification_model {castor} | R Documentation |
Simulate a deterministic uniform speciation/extinction model.
Description
Simulate a speciation/extinction cladogenic model for diversity over time, in the derministic limit. Speciation (birth) and extinction (death) rates can each be constant or power-law functions of the number of extant species. For example,
where is the birth rate,
is the intercept,
is the power-law factor,
is the current number of extant species and
is the power-law exponent. Optionally, the model can account for incomplete taxon sampling (rarefaction of tips) and for the effects of collapsing a tree at a non-zero resolution (i.e. clustering closely related tips into a single tip).
Usage
simulate_diversification_model( times,
parameters = list(),
added_rates_times = NULL,
added_birth_rates_pc = NULL,
added_death_rates_pc = NULL,
added_periodic = FALSE,
start_time = NULL,
final_time = NULL,
start_diversity = 1,
final_diversity = NULL,
reverse = FALSE,
include_coalescent = FALSE,
include_event_rates = FALSE,
include_Nevents = FALSE,
max_runtime = NULL)
Arguments
times |
Numeric vector, listing the times for which to calculate diversities, as predicted by the model. Values must be in ascending order. |
parameters |
A named list specifying the birth-death model parameters, with one or more of the following entries:
|
added_rates_times |
Numeric vector, listing time points (in ascending order) for a custom per-capita birth and/or death rates time series (see |
added_birth_rates_pc |
Numeric vector of the same size as |
added_death_rates_pc |
Numeric vector of the same size as |
added_periodic |
Logical, indicating whether |
start_time |
Numeric. Start time of the tree (<= |
final_time |
Numeric. Final (ending) time of the tree (>= |
start_diversity |
Numeric. Total diversity at |
final_diversity |
Numeric. Total diversity at |
reverse |
Logical. If |
include_coalescent |
Logical, specifying whether the diversity corresponding to a coalescent tree (i.e. the tree spanning only extant tips) should also be calculated. If |
include_event_rates |
Logical. If |
include_Nevents |
Logical. If |
max_runtime |
Numeric. Maximum runtime (in seconds) allowed for the simulation. If this time is surpassed, the simulation aborts. |
Details
The simulation is deterministic, meaning that diversification is modeled using ordinary differential equations, not as a stochastic process. The simulation essentially computes the deterministic diversity over time, not an actual tree. For stochastic cladogenic simulations yielding a random tree, see generate_random_tree
and simulate_dsse
.
In the special case where per-capita birth and death rates are constant (i.e. and
for birth and death rates), this function uses an explicit analytical solution to the underlying differential equations, and is thus much faster than in the general case.
If rarefaction<1
and resolution>0
, collapsing of closely related tips (at the resolution specified) is assumed to take place prior to rarefaction (i.e., subsampling applies to the already collapsed tips).
Value
A named list with the following elements:
success |
Logical, indicating whether the simulation was successful. If the simulation aborted due to runtime constraints (option |
total_diversities |
Numeric vector of the same size as |
coalescent_diversities |
Numeric vector of the same size as |
birth_rates |
Numeric vector of the same size as |
death_rates |
Numeric vector of the same size as |
Nbirths |
Numeric vector of the same size as |
Ndeaths |
Numeric vector of the same size as |
Author(s)
Stilianos Louca
See Also
generate_random_tree
,
count_lineages_through_time
Examples
# Generate a tree
max_time = 100
parameters = list(birth_rate_intercept = 10,
birth_rate_factor = 0,
birth_rate_exponent = 0,
death_rate_intercept = 0,
death_rate_factor = 0,
death_rate_exponent = 0,
resolution = 20,
rarefaction = 0.5)
generator = generate_random_tree(parameters,max_time=max_time)
tree = generator$tree
final_total_diversity = length(tree$tip.label)+generator$Nrarefied+generator$Ncollapsed
# Calculate diversity-vs-time curve for the tree
times = seq(from=0,to=0.99*max_time,length.out=50)
tree_diversities = count_lineages_through_time(tree, times=times)$lineages
# simulate diversity curve based on deterministic model
simulation = simulate_diversification_model(times,
parameters,
reverse=TRUE,
final_diversity=final_total_diversity,
include_coalescent=TRUE)
model_diversities = simulation$coalescent_diversities
# compare diversities in the tree to the simulated ones
plot(tree_diversities,model_diversities,xlab="tree diversities",ylab="simulated diversities")
abline(a=0,b=1,col="#A0A0A0") # show diagonal for reference