sample.species: Constant and time-dependent rate species sampling

Description Usage Arguments Value Author(s) Examples

View source: R/sample.species.R

Description

Generates a vector of occurrence times for a species in a simulation using a a Poisson process. Allows for the Poisson rate to be (1) a constant or (2) a function of time. For sampling of more than one species and/or taking into account species age instead of absolute time, see sample.clade and sample.adpp. sample.clade also allows for more flexibility options, see make.rate. Note that while the Poisson process occurs in forward time, we return (both in birth-death functions and here) results in backwards time, so that time is inverted using tMax both at the beginning and end of sample.species.

Usage

1
sample.species(sim, rr, tMax, S)

Arguments

sim

A sim object, usually an output of bd.sim.

rr

Sampling rate (per species per million years) over time. It can be a numeric describing a constant rate or a function(t) describing the variation in sampling over time. For more flexibility on sampling, see make.rate for creating more complex rates. Note that rr should always be greater than or equal to zero.

tMax

The maximum simulation time, used by rexp.var. A sampling time greater than tMax would mean the occurrence is sampled after the present, so for consistency we require this argument. This is also required to ensure time follows the correct direction both in the Poisson process and in the return.

S

The species number to be sampled. Since sample.species will be called by a wrapper using lapply, it is through S that we apply this function.

Value

A vector of occurrence times for that species.

Author(s)

Bruno do Rosario Petrucci and Matheus Januario.

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
###
# let us start with a linear increase in preservation rate

# simulate a group
sim <- bd.sim(n0 = 1, pp = 0.1, qq = 0.1, tMax = 10)

# in case first simulation was short-lived
while ((sim$TS[1] - ifelse(is.na(sim$TE[1]), 0, sim$TE[1])) < 10) {
  sim <- bd.sim(n0 = 1, pp = 0.1, qq = 0.1, tMax = 10)
}

# we will need to get exact durations for some examples, so
sim$TE[sim$EXTANT] = 0
# this is necessary since the default is to have NA for extant species

# preservation function
r <- function(t) {
  return(1 + 0.25*t)
}

# time
t <- seq(0, 10, by = 0.1)

# visualizing from the past to the present
plot(x = t, y = rev(r(t)), main="Simulated preservation", type = "l",
     xlab = "Mya", ylab = "preservation rate",
     xlim = c(10, sim$TE[1]))

# sample
occs <- sample.species(sim = sim, rr = r, tMax = 10, S = 1)

# check histogram
hist(occs,
     xlim = c(10, sim$TE[1]),
     xlab = "Mya")
lines(t, rev(r(t)))

###
# now let us try a step function

# simulate a group
sim <- bd.sim(n0 = 1, pp = 0.1, qq = 0.1, tMax = 10)

# in case first simulation was short lived
while ((sim$TS[1] - ifelse(is.na(sim$TE[1]), 0, sim$TE[1])) < 10) {
  sim <- bd.sim(n0 = 1, pp = 0.1, qq = 0.1, tMax = 10)
}

# we will need to get exact durations for some examples, so
sim$TE[sim$EXTANT] = 0
# this is necessary since the default is to have NA for extant species

# we can create the sampling rate here from a few vectors

# rates
rList <- c(1, 3, 0.5)

# rate shift times -  this could be c(10, 6, 2)
# and would produce the same function
rShifts <- c(0, 4, 8)

# create the rate to visualize it
r <- make.rate(rList, tMax = 10, fShifts = rShifts)

# time
t <- seq(0, 10, by = 0.1)

# visualizing the plot from past to present
plot(x = t, y = rev(r(t)), main = "Simulated preservation", type = "l",
     xlab = "Mya", ylab = "preservation rate",
     xlim = c(10, sim$TE[1]))

# sample
occs <- sample.species(sim = sim, rr = r, tMax = 10, S = 1)

# check histogram
hist(occs,
     xlim = c(10, sim$TE[1]),
     xlab = "Mya")

# frontiers of each regime
abline(v = 10 - rShifts, col = "red")

###
# we can create a step function in a different way as well

# simulate a group
sim <- bd.sim(n0 = 1, pp = 0.1, qq = 0.1, tMax = 10)

# in case first simulation was short-lived
while ((sim$TS[1] - ifelse(is.na(sim$TE[1]), 0, sim$TE[1])) < 10) {
  sim <- bd.sim(n0 = 1, pp = 0.1, qq = 0.1, tMax = 10)
}

# we will need to get exact durations for some examples, so
sim$TE[sim$EXTANT] = 0
# this is necessary since the default is to have NA for extant species

# preservation function
r <- function(t) {
  ifelse(t < 4, 1,
         ifelse(t < 8, 3, 0.5))
}
# note how this function should be exactly the same as the previous one

# time
t <- seq(0, 10, by = 0.1)

# visualizing the plot from past to present
plot(x = t, y = rev(r(t)), main = "Simulated preservation", type = "l",
     xlab = "Mya", ylab = "preservation rate",
     xlim = c(10, sim$TE[1]))

# sample
occs <- sample.species(sim = sim, rr = r, tMax = 10, S = 1)

# check histogram
hist(occs,
     xlim = c(10, sim$TE[1]),
     xlab = "Mya")
abline(v = 10 - rShifts, col = "red")

###
# finally we could generate sampling dependent on temperature

# simulate a group
sim <- bd.sim(n0 = 1, pp = 0.1, qq = 0.1, tMax = 10)

# in case first simulation was short-lived
while ((sim$TS[1] - ifelse(is.na(sim$TE[1]), 0, sim$TE[1])) < 10) {
  sim <- bd.sim(n0 = 1, pp = 0.1, qq = 0.1, tMax = 10)
}

# we will need to get exact durations for some examples, so
sim$TE[sim$EXTANT] = 0
# this is necessary since the default is to have NA for extant species

# preservation function dependent on temperature
r_t <- function(t, env) {
  return(0.25*env)
}

# get the temperature data
data(temp)

# final preservation
r <- make.rate(r_t, envF = temp)

# visualizing the plot from past to present
plot(x = t, y = rev(r(t)), main = "Simulated preservation", type = "l",
     xlab = "Mya", ylab = "preservation rate",
     xlim = c(10, ifelse(is.na(sim$TE[1]), 0, sim$TE[1])))

# sample
occs <- sample.species(sim = sim, rr = r, tMax = 10, S = 1)

# check histogram
hist(occs,
     xlim = c(10, ifelse(is.na(sim$TE[1]), 0, sim$TE[1])),
     xlab = "Mya")
lines(t, rev(r(t)))

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