Generate a random timetree via simulation of a Poissonian speciation/extinction (birth/death) process. New species are added (born) by splitting of a randomly chosen extant tip. The tree-wide birth and death rates of tips can each be constant or power-law functions of the number of extant tips. For example, $$ B = I + F\cdot N^E, $$ where \(B\) is the tree-wide birth rate (species generation rate), \(I\) is the intercept, \(F\) is the power-law factor, \(N\) is the current number of extant tips and \(E\) is the power-law exponent. Optionally, the per-capita (tip-specific) birth and death rates can be extended by adding a custom time series provided by the user.
generate_random_tree(parameters = list(),
max_tips = NULL,
max_time = NULL,
max_time_eq = NULL,
coalescent = TRUE,
as_generations = FALSE,
Nsplits = 2,
added_rates_times = NULL,
added_birth_rates_pc = NULL,
added_death_rates_pc = NULL,
added_periodic = FALSE,
tip_basename = "",
node_basename = NULL,
edge_basename = NULL,
include_birth_times = FALSE,
include_death_times = FALSE)A named list specifying the birth-death model parameters, with one or more of the following entries:
birth_rate_intercept: Non-negative number. The intercept of the Poissonian rate at which new species (tips) are added. In units 1/time. By default this is 0.
birth_rate_factor:
Non-negative number. The power-law factor of the Poissonian rate at which new species (tips) are added. In units 1/time. By default this is 0.
birth_rate_exponent:
Numeric. The power-law exponent of the Poissonian rate at which new species (tips) are added. Unitless. By default this is 1.
death_rate_intercept:
Non-negative number. The intercept of the Poissonian rate at which extant species (tips) go extinct. In units 1/time. By default this is 0.
death_rate_factor:
Non-negative number. The power-law factor of the Poissonian rate at which extant species (tips) go extinct. In units 1/time. By default this is 0.
death_rate_exponent:
Numeric. The power-law exponent of the Poissonian rate at which extant species (tips) go extinct. Unitless. By default this is 1.
resolution:
Non-negative numeric, specifying the resolution (in time units) at which to collapse the final tree by combining closely related tips. Any node whose age is smaller than this threshold, will be represented by a single tip. Set resolution=0 to not collapse tips (default).
rarefaction:
Numeric between 0 and 1. Rarefaction to be applied to the final tree (fraction of random tips kept in the tree). Note that if coalescent==FALSE, rarefaction may remove both extant as well as extinct clades. Set rarefaction=1 to not perform any rarefaction (default).
Maximum number of tips of the tree to be generated. If coalescent=TRUE, this refers to the number of extant tips. Otherwise, it refers to the number of extinct + extant tips. If NULL or <=0, the number of tips is unlimited (so be careful).
Maximum duration of the simulation. If NULL or <=0, this constraint is ignored.
Maximum duration of the simulation, counting from the first point at which speciation/extinction equilibrium is reached, i.e. when (birth rate - death rate) changed sign for the first time. If NULL or <0, this constraint is ignored.
Logical, specifying whether only the coalescent tree (i.e. the tree spanning the extant tips) should be returned. If coalescent==FALSE and the death rate is non-zero, then the tree may include non-extant tips (i.e. tips whose distance from the root is less than the total time of evolution). In that case, the tree will not be ultrametric.
Logical, specifying whether edge lengths should correspond to generations. If FALSE, then edge lengths correspond to time.
Integer greater than 1. Number of child-tips to generate at each diversification event. If set to 2, the generated tree will be bifurcating. If >2, the tree will be multifurcating.
Numeric vector, listing time points (in ascending order) for the custom per-capita birth and/or death rates time series (see added_birth_rates_pc and added_death_rates_pc below). Can also be NULL, in which case the custom time series are ignored.
Numeric vector of the same size as added_rates_times, listing per-capita birth rates to be added to the power law part. Can also be NULL, in which case this option is ignored and birth rates are purely described by the power law.
Numeric vector of the same size as added_rates_times, listing per-capita death rates to be added to the power law part. Can also be NULL, in which case this option is ignored and death rates are purely described by the power law.
Logical, indicating whether added_birth_rates_pc and added_death_rates_pc should be extended periodically if needed (i.e. if not defined for the entire simulation time). If FALSE, added birth & death rates are extended with zeros.
Character. Prefix to be used for tip labels (e.g. "tip."). If empty (""), then tip labels will be integers "1", "2" and so on.
Character. Prefix to be used for node labels (e.g. "node."). If NULL, no node labels will be included in the tree.
Character. Prefix to be used for edge labels (e.g. "edge."). Edge labels (if included) are stored in the character vector edge.label. If NULL, no edge labels will be included in the tree.
Logical. If TRUE, then the times of speciation events (in order of occurrence) will also be returned.
Logical. If TRUE, then the times of extinction events (in order of occurrence) will also be returned.
A named list with the following elements:
Logical, indicating whether the tree was successfully generated. If FALSE, the only other value returned is error.
A rooted bifurcating (if Nsplits==2) or multifurcating (if Nsplits>2) tree of class "phylo", generated according to the specified birth/death model. If coalescent==TRUE or if all death rates are zero, and only if as_generations==FALSE, then the tree will be ultrametric. If as_generations==TRUE and coalescent==FALSE, all edges will have unit length.
Numeric, giving the time at which the tree's root was first split during the simulation. Note that if coalescent==TRUE, this may be later than the first speciation event during the simulation.
Numeric, giving the final time at the end of the simulation. Note that if coalescent==TRUE, then this may be greater than the total time span of the tree (since the root of the coalescent tree need not correspond to the first speciation event).
Numeric, giving the first time where the sign of (death rate - birth rate) changed from the beginning of the simulation, i.e. when speciation/extinction equilibrium was reached. May be infinite if the simulation stoped before reaching this point.
Integer vector, listing indices of extant tips in the tree. If coalescent==TRUE, all tips will be extant.
Total number of birth events (speciations) that occurred during tree growth. This may be lower than the total number of tips in the tree if death rates were non-zero and coalescent==TRUE, or if Nsplits>2.
Total number of deaths (extinctions) that occurred during tree growth.
Number of tips removed from the tree while collapsing at the resolution specified.
Number of tips removed from the tree due to rarefaction.
Numeric vector, listing the times of speciation events during tree growth, in order of occurrence. Note that if coalescent==TRUE, then speciation_times may be greater than the phylogenetic distance to the coalescent root.
Numeric vector, listing the times of extinction events during tree growth, in order of occurrence. Note that if coalescent==TRUE, then speciation_times may be greater than the phylogenetic distance to the coalescent root.
Character, containing an explanation of ther error that occurred. Only included if success==FALSE.
If max_time==NULL, then the returned tree will always contain max_tips tips. In particular, if at any moment during the simulation the tree only includes a single extant tip, the death rate is temporarily set to zero to prevent the complete extinction of the tree. If max_tips==NULL, then the simulation is ran as long as specified by max_time. If neither max_time nor max_tips is NULL, then the simulation halts as soon as the time exceeds max_time or the number of tips (extant tips if coalescent is TRUE) exceeds max_tips. If max_tips!=NULL and Nsplits>2, then the last diversification even may generate fewer than Nsplits children, in order to keep the total number of tips within the specified limit.
If rarefaction<1 and resolution>0, collapsing of closely related tips (at the resolution specified) takes place prior to rarefaction (i.e., subsampling applies to the already collapsed tips).
Both the per-capita birth and death rates can be made into completely arbitrary functions of time, by setting all power-law coefficients to zero and providing custom time series added_birth_rates_pc and added_death_rates_pc.
D. J. Aldous (2001). Stochastic models and descriptive statistics for phylogenetic trees, from Yule to today. Statistical Science. 16:23-34.
M. Steel and A. McKenzie (2001). Properties of phylogenetic trees generated by Yule-type speciation models. Mathematical Biosciences. 170:91-112.
# NOT RUN {
# Simple speciation model
parameters = list(birth_rate_intercept=1)
tree = generate_random_tree(parameters,max_tips=100)$tree
# Exponential growth rate model
parameters = list(birth_rate_factor=1)
tree = generate_random_tree(parameters,max_tips=100)$tree
# }
Run the code above in your browser using DataLab