generate_tree_with_evolving_rates {castor} | R Documentation |

Generate a random phylogenetic tree via simulation of a Poissonian speciation/extinction (birth/death) process. New species are added (born) by splitting of a randomly chosen extant tip. Per-capita birth and death rates (aka. speciation and extinction rates) evolve under some stochastic process (e.g. Brownian motion) along each edge. Thus, the probability rate of a tip splitting or going extinct depends on the tip, with closely related tips having more similar per-capita birth and death rates.

generate_tree_with_evolving_rates(parameters = list(), rate_model = 'BM', max_tips = NULL, max_time = NULL, max_time_eq = NULL, coalescent = TRUE, as_generations = FALSE, tip_basename = "", node_basename = NULL, include_event_times = FALSE, include_rates = FALSE)

`parameters` |
A named list specifying the model parameters for the evolving birth/death rates. The precise entries expected depend on the chosen |

`rate_model` |
Character, specifying the model for the evolving per-capita birth/death rates. Must be one of the following: 'BM' (Brownian motion constrained to a finite interval via reflection), 'Mk' (discrete-state continuous-time Markov chain with fixed transition rates). |

`max_tips` |
Maximum number of tips of the tree to be generated. If |

`max_time` |
Maximum duration of the simulation. If |

`max_time_eq` |
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 |

`coalescent` |
Logical, specifying whether only the coalescent tree (i.e. the tree spanning the extant tips) should be returned. If |

`as_generations` |
Logical, specifying whether edge lengths should correspond to generations. If FALSE, then edge lengths correspond to time. |

`tip_basename` |
Character. Prefix to be used for tip labels (e.g. "tip."). If empty (""), then tip labels will be integers "1", "2" and so on. |

`node_basename` |
Character. Prefix to be used for node labels (e.g. "node."). If |

`include_event_times` |
Logical. If |

`include_rates` |
Logical. If |

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 `rate_model=='BM'`

, then per-capita birth rates (speciation rates) and per-capita death rates (extinction rates) evolve according to Brownian Motion, constrained to a finite interval via reflection. Note that speciation and extinction rates are only updated at branching points, i.e. during speciation events, while waiting times until speciation/extinction are based on rates at the previous branching point. The argument `parameters`

should be a named list including one or more of the following elements:

`birth_rate_diffusivity`

: Non-negative number. Diffusivity constant for the Brownian motion model of the evolving per-capita birth rate. In units 1/time^3. See`simulate_bm_model`

for an explanation of the diffusivity parameter.`min_birth_rate_pc`

: Non-negative number. The minimum allowed per-capita birth rate of a clade. In units 1/time. By default this is 0.`max_birth_rate_pc`

: Non-negative number. The maximum allowed per-capita birth rate of a clade. In units 1/time. By default this is 1.`death_rate_diffusivity`

: Non-negative number. Diffusivity constant for the Brownian motion model of the evolving per-capita death rate. In units 1/time^3. See`simulate_bm_model`

for an explanation of the diffusivity parameter.`min_death_rate_pc`

: Non-negative number. The minimum allowed per-capita death rate of a clade. In units 1/time. By default this is 0.`max_death_rate_pc`

: Non-negative number. The maximum allowed per-capita death rate of a clade. In units 1/time. By default this is 1.`root_birth_rate_pc`

: Non-negative number, between`min_birth_rate_pc`

and`max_birth_rate_pc`

, specifying the initial per-capita birth rate of the root. If left unspecified, this will be chosen randomly and uniformly within the allowed interval.`root_death_rate_pc`

: Non-negative number, between`min_death_rate_pc`

and`max_death_rate_pc`

, specifying the initial per-capita death rate of the root. If left unspecified, this will be chosen randomly and uniformly within the allowed interval.`rarefaction`

: Numeric between 0 and 1. Rarefaction to be applied at the end of the simulation (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.

If `rate_model=='Mk'`

, then speciation/extinction rates are determined by a tip's current "state", which evolves according to a continuous-time discrete-state Markov chain (Mk model) with constant transition rates. The argument `parameters`

should be a named list including one or more of the following elements:

`Nstates`

: Number of possible discrete states a tip can have. For example, if`Nstates`

then this corresponds to the common Binary State Speciation and Extinction (BiSSE) model (Maddison et al., 2007). By default this is 1.`state_birth_rates`

: Numeric vector of size Nstates, listing the per-capita birth rate (speciation rate) at each state. Can also be a single number (all states have the same birth rate).`state_death_rates`

: Numeric vector of size Nstates, listing the per-capita death rate (extinction rate) at each state. Can also be a single number (all states have the same death rate).`transition_matrix`

: 2D numeric matrix of size Nstates x Nstates. Transition rate matrix for the Markov chain model of birth/death rate evolution.`start_state`

: Integer within 1,..,`Nstates`

, specifying the initial state of the first created lineage. If left unspecified, this is chosen randomly and uniformly among all possible states.`rarefaction`

: Same as when`rate_model=='BM'`

.

Note: The option `rate_model=='Mk'`

is deprecated and included for backward compatibility purposes only. To generate a tree with Markov transitions between states (known as Multiple State Speciation and Extinction model), use the command `simulate_dsse`

instead.

A named list with the following elements:

`success` |
Logical, indicating whether the simulation was successful. If |

`tree` |
A rooted bifurcating tree of class "phylo", generated according to the specified birth/death model. If |

`root_time` |
Numeric, giving the time at which the tree's root was first split during the simulation.
Note that if |

`final_time` |
Numeric, giving the final time at the end of the simulation. If |

`equilibrium_time` |
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. |

`Nbirths` |
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 |

`Ndeaths` |
Total number of deaths (extinctions) that occurred during tree growth. |

`birth_times` |
Numeric vector, listing the times of speciation events during tree growth, in order of occurrence. Note that if |

`death_times` |
Numeric vector, listing the times of extinction events during tree growth, in order of occurrence. Note that if |

`birth_rates_pc` |
Numeric vector, listing the per-capita birth rate of each tip and node in the tree. The length of an edge in the tree was thus drawn from an exponential distribution with rate equal to the per-capita birth rate of the child tip or node. |

`death_rates_pc` |
Numeric vector, listing the per-capita death rate of each tip and node in the tree. |

`states` |
Integer vector of size Ntips+Nnodes, listing the discrete state of each tip and node in the tree. Only included if |

`start_state` |
Integer, specifying the initial state of the first created lineage (either provided during the function call, or generated randomly). Only included if |

`root_birth_rate_pc` |
Numeric, specifying the initial per-capita birth rate of the root (either provided during the function call, or generated randomly). Only included if |

`root_death_rate_pc` |
Numeric, specifying the initial per-capita death rate of the root (either provided during the function call, or generated randomly). Only included if |

Stilianos Louca

D. J. Aldous (2001). Stochastic models and descriptive statistics for phylogenetic trees, from Yule to today. Statistical Science. 16:23-34.

W. P. Maddison, P. E. Midford, S. P. Otto (2007). Estimating a binary character's effect on speciation and extinction. Systematic Biology. 56:701-710.

# Example 1 # Generate tree, with rates evolving under Brownian motion parameters = list(birth_rate_diffusivity = 1, min_birth_rate_pc = 1, max_birth_rate_pc = 2, death_rate_diffusivity = 0.5, min_death_rate_pc = 0, max_death_rate_pc = 1) simulation = generate_tree_with_evolving_rates(parameters, rate_model='BM', max_tips=1000, include_rates=TRUE) tree = simulation$tree Ntips = length(tree$tip.label) # plot per-capita birth & death rates of tips plot( x=simulation$birth_rates_pc[1:Ntips], y=simulation$death_rates_pc[1:Ntips], type='p', xlab="pc birth rate", ylab="pc death rate", main="Per-capita birth & death rates across tips (BM model)", las=1) ###################### # Example 2 # Generate tree, with rates evolving under a binary-state model Q = get_random_mk_transition_matrix(Nstates=2, rate_model="ER", max_rate=0.1) parameters = list(Nstates = 2, state_birth_rates = c(1,1.5), state_death_rates = 0.5, transition_matrix = Q) simulation = generate_tree_with_evolving_rates(parameters, rate_model='Mk', max_tips=1000, include_rates=TRUE) tree = simulation$tree Ntips = length(tree$tip.label) # plot distribution of per-capita birth rates of tips rates = simulation$birth_rates_pc[1:Ntips] barplot(table(rates)/length(rates), xlab="rate", main="Distribution of pc birth rates across tips (Mk model)")

[Package *castor* version 1.6.8 Index]