# var.rate.div: Expected diversity for general exponential rates In brpetrucci/paleobuddy: Simulating Diversification Dynamics

## 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 numberFor constant birth-death rates. A function of timeFor rates that vary with time. Note that this can be any function of time. A function of time and an environmental variableFor 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 vectorTo 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 < fShifts`), 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.