Description Usage Arguments Value Author(s) Examples

View source: R/sample.species.R

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`

.

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

`sim` |
A |

`rr` |
Sampling rate (per species per million years) over time. It can be
a |

`tMax` |
The maximum simulation time, used by |

`S` |
The species number to be sampled. Since |

A vector of occurrence times for that species.

Bruno do Rosario Petrucci and Matheus Januario.

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)))
``` |

Embedding an R snippet on your website

Add the following code to your website.

For more information on customizing the embed code, read Embedding Snippets.