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

Description Usage Arguments Value Examples

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.

1	var.rate.div(rate, t, n0 = 1, tMax = NULL, envRate = NULL, rateShifts = NULL)

`rate`	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 `envRate` will result in an error. A numeric vector To create step function rates. Note this must be accompanied by a corresponding vector of rate shift times, `rateShifts`.
`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 species given a rate `rate` and an initial number of parents `n0`, so in a biological context `rate` 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 `rateShifts` runs the correct way.
`envRate`	A `data.frame` representing a time-series, usually an environmental variable (e.g. CO2, temperature, etc) varying 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 the user supplies `envRate` without `rate` being a function of two variables. `paleobuddy` has two environmental data frames, `temp` and `co2`. One can check `RPANDA` for more examples. Note that, since simulation functions are run in forward-time (i.e. with `0` being the origin time of the simulation), the time column of `envRate` is assumed to do so as well, so that the row corresponding to `t = 0` is assumed to be the value of the time-series when the simulation starts, and `t = tMax` is assumed to be its value when the simulation ends (the present). Acknowledgements: The strategy to transform a function of `t` and `envRate` into a function of `t` only using `envRate` was adapted from `RPANDA`.
`rateShifts`	A vector indicating the time of rate shifts in a step function. The first element must be the first or last time point for the simulation, i.e. `0` or `tMax`. Since functions in paleobuddy run from `0` to `tMax`, if `rateShifts` runs from past to present (meaning `rateShifts[2] < rateShifts[1]`), we take `tMax - rateShifts` as the shifts vector. Note that supplying `rateShifts` when `rate` is not a numeric vector of the same length will result in an error.

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

# let us first create a vector of times to use in these examples
time <- seq(0, 50, 0.1)

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

# make the rate
r <- make.rate(0.5)

# plot it
plot(time, rep(r, length(time)), ylab = "Diversification rate",
     xlab = "Time (Mya)", xlim = c(50, 0), type = 'l')

# get expected diversity
div <- var.rate.div(rate, time)

# plot it
plot(time, rev(div), ylab = "Expected number of species",
     xlab = "Time (Mya)", xlim = c(50, 0), type = 'l')

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

# make the rate
r <- make.rate(rate)

# plot it
plot(time, rev(r(time)), ylab = "Diversification rate",
     xlab = "Time (Mya)", xlim = c(50, 0), type = 'l')
# negative values are ok since this represents a diversification rate

# get expected diversity
div <- var.rate.div(rate, time)

# plot it
plot(time, rev(div), ylab = "Expected number of species",
     xlab = "Time (Mya)", xlim = c(50, 0), type = 'l')

###
# remember: rate is the diversification rate!

# we can create speciation...
lambda <- function(t) {
  return(0.5 - 0.01*t)
}

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

# ...and get rate as diversification
rate <- function(t) {
  return(lambda(t) - mu(t))
}

# make the rate
r <- make.rate(rate)

# plot it
plot(time, rev(r(time)), ylab = "Diversification rate",
     xlab = "Time (Mya)", xlim = c(50, 0), type = 'l')

# get expected diversity
div <- var.rate.div(rate, time)

# plot it
plot(time, rev(div), ylab = "Expected number of species",
     xlab = "Time (Mya)", xlim = c(50, 0), type = 'l')

###
# we can use ifelse() to make a step function like this
rate <- function(t) {
  return(ifelse(t < 2, 0.2,
          ifelse(t < 3, 0.4,
           ifelse(t < 5, -0.2, 0.5))))
}

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

# make the rate
r <- make.rate(rate)

# plot it
plot(time, rev(r(time)), ylab = "Diversification rate",
     xlab = "Time (Mya)", xlim = c(10, 0), type = 'l')
# negative rates is ok since this represents a diversification rate

# get expected diversity
div <- var.rate.div(rate, time)

# plot it
plot(time, rev(div), ylab = "Expected number of species",
     xlab = "Time (Mya)", xlim = c(10, 0), type = 'l')

# 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 a shifts vector into a function of time

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

# rates vector
rateList <- c(0.2, 0.4, -0.2, 0.5)

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

# make the rate
r <- make.rate(rateList, tMax = 10, rateShifts = rateShifts)

# plot it
plot(time, rev(r(time)), ylab = "Diversification rate",
     xlab = "Time (Mya)", xlim = c(10, 0), type = 'l')
# negative rates is ok since this represents a diversification rate

# get expected diversity
div <- var.rate.div(rateList, time, tMax = 10, rateShifts = rateShifts)

# plot it
plot(time, rev(div), ylab = "Expected number of species",
     xlab = "Time (Mya)", xlim = c(10, 0), type = 'l')
 
###
# finally let us see what we can do with environmental variables

# get the temperature data
data(temp)

# diversification
rate <- function(t, env) {
  return(0.2 + 2*exp(-0.25*env))
}

# make the rate
r <- make.rate(rate, tMax = tMax, envRate = temp)

# plot it
plot(time, rev(r(time)), ylab = "Diversification rate",
     xlab = "Time (Mya)", xlim = c(10, 0), type = 'l')

# get expected diversity
div <- var.rate.div(rate, time, tMax = tMax, envRate = temp)

# plot it
plot(time, rev(div), ylab = "Expected number of species",
     xlab = "Time (Mya)", xlim = c(10, 0), type = 'l')

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

# diversification
rate <- function(t, env) {
  return(0.02 * env - 0.01 * t)
}

# make the rate
r <- make.rate(rate, tMax = tMax, envRate = temp)

# plot it
plot(time, rev(r(time)), ylab = "Diversification rate",
     xlab = "Time (Mya)", xlim = c(10, 0), type = 'l')

# get expected diversity
div <- var.rate.div(rate, time, tMax = tMax, envRate = temp)

# plot it
plot(time, rev(div), ylab = "Expected number of species",
     xlab = "Time (Mya)", xlim = c(10, 0), type = 'l')
  
###
# as mentioned above, we could also use ifelse() to construct a step 
# function that is modulated by temperature

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

# make the rate
r <- make.rate(rate, tMax = tMax, envRate = temp)

# plot it
plot(time, rev(r(time)), ylab = "Diversification rate",
     xlab = "Time (Mya)", xlim = c(10, 0), type = 'l')


# get expected diversity
div <- var.rate.div(rate, time, tMax = tMax, envRate = temp)

# plot it
plot(time, rev(div), ylab = "Expected number of species",
     xlab = "Time (Mya)", xlim = c(10, 0), type = 'l')


# takes a bit long so we set it to not run, but the user
# should feel free to explore this and other scenarios