simulate_diversification_model {castor}  R Documentation 
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 powerlaw functions of the number of extant species. For example,
B = I + F\cdot N^E,
where B is the birth rate, I is the intercept, F is the powerlaw factor, N is the current number of extant species and E is the powerlaw exponent. Optionally, the model can account for incomplete taxon sampling (rarefaction of tips) and for the effects of collapsing a tree at a nonzero resolution (i.e. clustering closely related tips into a single tip).
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)
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 birthdeath model parameters, with one or more of the following entries:

added_rates_times 
Numeric vector, listing time points (in ascending order) for a custom percapita 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. 
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 percapita birth and death rates are constant (i.e. I=0 and E=1 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).
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 
Stilianos Louca
generate_random_tree
,
count_lineages_through_time
# 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 diversityvstime 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