bd.sim.general: Non-constant rate Birth-Death simulation

Description Usage Arguments Value Author(s) Examples

View source: R/bd.sim.general.R

Description

Simulates a species birth-death process with general rates for any number of starting species. Allows for the speciation/extinction rate to be (1) a constant, or (2) a function of time. Allows for constraining results on the number of species at the end of the simulation, either total or extant. The function can also take an optional shape argument to generate age-dependence on speciation and/or extinction, assuming a Weibull distribution as a model of age-dependence. Returns an object containing vectors of speciation times, extinction times, parents (= species' mother species) and status at the end of the simulation (extant or not) for each species in the simulation. It may return true extinction times or simply information on whether species lived after the maximum simulation time. For constant rate simulations, see bd.sim.constant. For a function that unites all scenarios, see bd.sim. bd.sim also allows for extra inputs, creating a time-dependent only rate internally through make.rate. For similar flexibility, use make.rate to generate the desired rate. Please note while time runs from 0 to tMax in the simulation, it returns speciation/extinction times as tMax (origin of the group) to 0 (the "present" and end of simulation), so as to conform to other packages in the literature.

Usage

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
bd.sim.general(
  n0,
  pp,
  qq,
  tMax,
  pShape = NULL,
  qShape = NULL,
  nFinal = c(0, Inf),
  extOnly = FALSE,
  trueExt = FALSE
)

Arguments

n0

Initial number of species. Usually 1, in which case the simulation will describe the full diversification of a monophyletic lineage. Notice that when pp is less than or equal to qq, many simulations will go extinct before speciating even once. One way of generating simulations in this case is to increase n0. The user should notice this will simulate the diversification of a paraphyletic group.

pp

Function to hold the speciation rate over time. It will either be interpreted as an exponential rate, or a Weibull scale if pShape != NULL.

qq

Similar to above, but for the extinction rate.

Note: this function is meant to be called by bd.sim, so it neither allows for as much flexibility, nor calls make.rate. If the user wishes to use bd.sim.general with environmental or step-function rates, they can generate the rate with make.rate and supply it to the function.

tMax

Ending time of simulation, in million years after the clade origin. Any species still living after tMax is considered extant, and any species that would be generated after tMax is not born.

pShape

Shape of the age-dependency in speciation rate. This will be equal to the shape parameter in a Weibull distribution: when smaller than one, speciation rate will decrease along each species' age (negative age-dependency). When larger than one, speciation rate will increase along each species's age (positive age-dependency). Default is NULL, equivalent to an age-independent process. For pShape != NULL (including when equal to one), pp will be considered a scale (= 1/rate), and rexp.var will draw a Weibull distribution instead of an exponential. This means Weibull(rate, 1) = Exponential(1/rate). Note that even when pShape != NULL, pp may still be time-dependent.

qShape

Similar to pShape, but for the extinction rate.

Note: Time-varying shape is implemented, so one could have pShape or qShape be a function of time. It is not thoroughly tested, however, so it may be prudent to wait for a future release where this feature is well established.

nFinal

A vector of length 2, indicating an interval of acceptable number of species at the end of the simulation. Default value is c(0, Inf), so that any number of species (including zero, the extinction of the whole clade) is accepted. If different from default value, the process will run until the number of total species reaches a number in the interval nFinal. Please note that extOnly can modify the meaning of this parameter. If extOnly = TRUE, the function will run repeatedly until the number of extant species at the end of the simulation lies within the nFinal interval. Note that using values other than the default for nFinal might condition simulation results.

extOnly

A logical indicating whether nFinal should be taken as an interval of the number of total or extant species during the simulation. If TRUE, the function will run repeatedly until the number of extant species at the end of the simulation lies within the nFinal interval. If FALSE (as default), it will run until the total number of species generated lies within that interval.

Note: The function returns NA if it runs for more than 100000 iterations without fulfilling the requirements of nFinal and extOnly.

trueExt

A logical used for bd.sim.general, indicating whether it should return true or truncated extinction times. When TRUE, time of extinction of extant species will be the true time, otherwise it will be NA if a species is alive at the end of the simulation.

Value

A list of vectors, as follows

TE

List of extinction times, with 0 as the time of extinction for extant species.

TS

List of speciation times, with NA as the time of speciation for species that started the simulation.

PAR

List of parents. Species that started the simulation have NA, while species that were generated during the simulation have their parent's number. Species are numbered as they are born.

EXTANT

List of logicals representing whether each species is extant.

Author(s)

Bruno do Rosario Petrucci.

Examples

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
# we can test a couple scenarios

###
# first, even though this is bd.sim.general, we can try constant rates

# initial number of species
n0 <- 1

# maximum simulation time
tMax <- 40

# speciation
p <- 0.11

# extinction
q <- 0.08

# run the simulation
sim <- bd.sim.general(n0, p, q, tMax, nFinal = c(2, Inf))

# we can plot the phylogeny to take a look
if (requireNamespace("ape", quietly = TRUE)) {
  phy <- make.phylo(sim)
  ape::plot.phylo(phy)
}

###
# we can complicate things further with a linear function as a rate
# bd.sim.general takes longer so we run examples for 1000 replicates instead

# initial number of species
n0 <- 1

# maximum simulation time
tMax <- 40

# speciation
p <- function(t) {
  return(0.03 + 0.005*t)
}

# extinction
q <- 0.05

# run the simulation
sim <- bd.sim.general(n0, p, q, tMax, nFinal = c(2, Inf))

# we can plot the phylogeny to take a look
if (requireNamespace("ape", quietly = TRUE)) {
  phy <- make.phylo(sim)
  ape::plot.phylo(phy)
}

###
# we can also create a step function

# initial number of species
n0 <- 1

# maximum simulation time
tMax <- 40

# speciation rate
p <- function(t) {
  return(0.03 + 0.005*t)
}

# vector of extinction rates
qList <- c(0.06, 0.09, 0.11)

# vector of shift times. Note qShifts could be c(40, 20, 10) for
# identical results
qShifts <- c(0, 15, 25)

# let us take a look at how make.rate will make it a step function
q <- make.rate(qList, tMax = tMax, fShifts = qShifts)

# and plot it
plot(seq(0, tMax, 0.1), q(seq(0, tMax, 0.1)), type = 'l',
     main = "Extintion rate as a step function", xlab = "Time (My)",
     ylab = "Rate (species/My)")

# a different way to define the same extinction function
q <- function(t) {
  ifelse(t < 15, 0.06,
         ifelse(t < 25, 0.09, 0.11))
}

# run the simulation
sim <- bd.sim.general(n0, p, q, tMax, nFinal = c(2, Inf))
# we could instead have used q made with make.rate
# that is, however, much slower

# we can plot the phylogeny to take a look
if (requireNamespace("ape", quietly = TRUE)) {
  phy <- make.phylo(sim)
  ape::plot.phylo(phy)
}

###
# another feature to add is age dependency
n0 <- 1

# maximum simulation time
tMax <- 40

# speciation
p <- 0.1

# extinction - a Weibull scale
q <- 10

# extinction shape
qShape <- 1

# run simulations
sim <- bd.sim.general(n0, p, q, tMax, qShape = qShape, nFinal = c(2, Inf))

# we can plot the phylogeny to take a look
if (requireNamespace("ape", quietly = TRUE)) {
  phy <- make.phylo(sim)
  ape::plot.phylo(phy)
}

###
# scale can be time-dependent
n0 <- 1

# maximum simulation time
tMax <- 40

# speciation
p <- 0.15

# extinction - a Weibull scale
q <- function(t) {
  return(8 + 0.05*t)
}

# extinction shape
qShape <- 1

# run simulations
sim <- bd.sim.general(n0, p, q, tMax, qShape = qShape, nFinal = c(2, Inf))

# we can plot the phylogeny to take a look
if (requireNamespace("ape", quietly = TRUE)) {
  phy <- make.phylo(sim)
  ape::plot.phylo(phy)
}

###
# finally, we could have environmental dependency on a rate
# initial number of species
n0 <- 1

# maximum simulation time
tMax <- 40

# temperature-dependent speciation
p_t <- function(t, temp) {
 return(0.025*exp(0.1*temp))
}

# extinction
q <- 0.075

# get the temperature data
data(temp)

# speciation
p <- make.rate(p_t, envF = temp)

# run simulations
sim <- bd.sim.general(n0, p, q, tMax, nFinal = c(2, Inf))

# we can plot the phylogeny to take a look
if (requireNamespace("ape", quietly = TRUE)) {
  phy <- make.phylo(sim)
  ape::plot.phylo(phy)
}

# note nFinal has to be sensible
## Not run: 
# this would return a warning, since it is virtually impossible to get 100
# species at a process with diversification rate -0.09 starting at n0 = 1
sim <- bd.sim.general(1, pp = 0.01, qq = 1, tMax = 100, nFinal = c(100, Inf))

## End(Not run)

brpetrucci/paleobuddy documentation built on Aug. 8, 2020, 2:03 a.m.