var.rate.div: Expected diversity for general exponential rates

Description Usage Arguments Value Examples

View source: R/var.rate.div.R

Description

Calculates the expected species diversity on an interval given a (possibly time dependent) exponential rate. Takes as the base rate (1) a constant, (2) a function of time, (3) a function of time interacting with an environmental variable, or (4) a vector of numbers describing rates as a step function. Requires information regarding the maximum simulation time, and allows for optional extra parameters to tweak the baseline rate. For more information on the creation of the final rate, see make.rate.

Usage

1
var.rate.div(ff, t, n0 = 1, tMax = NULL, envF = NULL, fShifts = NULL)

Arguments

ff

The baseline function with which to make the rate. It can be a

A number

For constant birth-death rates.

A function of time

For rates that vary with time. Note that this can be any function of time.

A function of time and an environmental variable

For rates varying with time and an environmental variable, such as temperature. Note that supplying a function on more than one variable without an accompanying envF will result in an error.

A numeric vector

To create step function rates. Note this must be accompanied by a corresponding vector of shifts fShifts.

t

A time vector over which to consider the distribution.

n0

The initial number of species is by default 1, but one can change to any nonnegative number.

Note: var.rate.div will find the expected number of daughters given a rate ff and an initial number of parents n0, so in a biological context ff is diversification rate, not speciation (unless extinction is 0).

tMax

Ending time of simulation, in million years after the clade's origin. Needed to ensure fShifts runs the correct way.

envF

A data.frame representing the variation of an environmental variable (e.g. CO2, temperature, available niches, etc) with time. The first column of this data.frame must be time, and the second column must be the values of the variable. The function will return an error if supplying envF without ff being a function of two variables. Note paleobuddy has two environmental data frames, temp and co2. One can check RPANDA for more examples.

Acknowledgements: The strategy to transform a function of t and env into a function of t only using envF was adapted from RPANDA (see below).

fShifts

A vector indicating the time placement of rate shifts in a step function. The first element must be the first time point for the simulation. This may be 0 or tMax. Since functions in paleobuddy run from 0 to tMax, if fShifts runs from past to present (fShifts[2] < fShifts[1]), we take tMax - fShifts as the shifts vector. Note that supplying fShifts when ff is not a numeric vector of the same length will result in an error.

Value

A vector of the expected number of species per time point supplied in t, which can then be used to plot vs. t.

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
# let us first create a vector of times to use in these examples.
t <- seq(0, 50, 0.1)

###
# we can start simple: create a constant rate
ff <- 0.1

# set this up so we see rates next to diversity
par(mfrow = c(1,2))

# see how the rate looks
r <- make.rate(0.5)
plot(t, rep(r, length(t)), type = 'l')

# get the diversity and plot it
div <- var.rate.div(ff, t)
plot(t, div, type = 'l')

###
# something a bit more complex: a linear rate
ff <- function(t) {
  return(0.01*t)
}

# visualize the rate
r <- make.rate(ff)
plot(t, r(t), type = 'l')

# get the diversity and plot it
div <- var.rate.div(ff, t = t)
plot(t, div, type = 'l')

###
# remember: ff is diversity!

# we can create speciation...
pp <- function(t) {
  return(-0.01*t + 0.2)
}

# ...and extinction...
qq <- function(t) {
  return(0.01*t)
}

# ...and code ff as diversification
ff <- function(t) {
  return(pp(t) - qq(t))
}

# visualize the rate
r <- make.rate(ff)
plot(t, r(t), type = 'l')

# get diversity and plot it
div <- var.rate.div(ff, t, n0 = 2)
plot(t, div, type = 'l')

###
# remember: ff can be any time-varying function!

# such as a sine
ff <- function(t) {
  return(sin(t)*0.5)
}

# visualize the rate
r <- make.rate(ff)
plot(t, r(t), type = 'l')

# we can have any number of starting species
div <- var.rate.div(ff, t, n0 = 2)
plot(t, div, type = 'l')

###
# we can use ifelse() to make a step function like this
ff <- function(t) {
  return(ifelse(t < 2, 0.1,
          ifelse(t < 3, 0.3,
           ifelse(t < 5, -0.2, 0.05))))
}

# change t so things are faster
t <- seq(0, 10, 0.1)

# visualize the rate
r <- make.rate(ff)
plot(t, r(t), type = 'l')

# get the diversity and plot it
div <- var.rate.div(ff, t)
plot(t, div, type = 'l')

# important note: this method of creating a step function might be annoying,
# but when running thousands of simulations it will provide a much faster
# integration than when using our method of transforming a rates and shifts
# vector into a function of time

###
# ...which we can do as follows

# rates vector
ff <- c(0.1, 0.3, -0.2, 0.05)

# rate shifts vector
fShifts <- c(0, 2, 3, 5)

# visualize the rate
r <- make.rate(ff, tMax = 10, fShifts = fShifts)
plot(t, r(t),type = 'l')

# get the diversity and plot it
div <- var.rate.div(ff, t, tMax = 10, fShifts = fShifts)
plot(t, div, type = 'l')

# note the delay in running var.rate.div using this method. integrating a step
# function created using the methods in make.rate() is slow, as explained in
# the make.rate documentation)

# it is also impractical to supply a rate and a shifts vector and
# have an environmental dependency, so in cases where one looks to run
# more than a couple dozen simulations, and when one is looking to have a
# step function modified by an environmental variable, consider using ifelse()

# finally let us see what we can do with environmental variables

# get the temperature data
data(temp)

# diversification
ff <- function(t, env) {
  return(0.002*env)
}

# visualize the rate
r <- make.rate(ff, envF = temp)
plot(t, r(t), type = 'l')

# get diversity and plot it
div <- var.rate.div(ff, t, envF = temp)
plot(t, div, type = 'l')

###
# we can also have a function that depends on both time AND temperature

# diversification
ff <- function(t, env) {
  return(0.02 * env - 0.001 * t)
}

# visualize the rate
r <- make.rate(ff, envF = temp)
plot(t, r(t), type = 'l')

# get diversity and plot it
div <- var.rate.div(ff, t, envF = temp)
plot(t, div, type = 'l')
  
###
# as mentioned above, we could also use ifelse() to construct a step function
# that is modulated by temperature

# diversification
ff <- function(t, env) {
  return(ifelse(t < 2, 0.1 + 0.01*env,
          ifelse(t < 5, 0.2 - 0.005*env,
           ifelse(t < 8, 0.1 + 0.005*env, 0.08))))
}

# visualize the rate
r <- make.rate(ff, envF = temp)
plot(t, r(t), type = 'l')

# get diversity and plot it
div <- var.rate.div(ff, t, envF = temp)
plot(t, div, type = 'l')

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