Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Dev #7

Merged
merged 5 commits into from
Jul 20, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
6 changes: 6 additions & 0 deletions R/StratPal-package.R
Original file line number Diff line number Diff line change
@@ -0,0 +1,6 @@
#' @keywords internal
"_PACKAGE"

## usethis namespace: start
## usethis namespace: end
NULL
14 changes: 9 additions & 5 deletions R/apply_niche_pref.R
Original file line number Diff line number Diff line change
Expand Up @@ -3,13 +3,13 @@ apply_niche_pref = function(x, niche_def, gc){
#'
#' @title apply niche preference
#'
#' @param x events, e.g. fossil occurrences
#' @param niche_def function, specifying the niche along a gradient. should return 0 when taxon is outside of niche, and 1 when fully inside niche. Values between 0 and 1 are interpreted as probabilities.
#' @param x events, e.g. times/ages of fossil occurrences or their stratigraphic position.
#' @param niche_def function, specifying the niche along a gradient. Should return 0 when taxon is outside of niche, and 1 when fully inside niche. Values between 0 and 1 are interpreted as probabilities.
#' @param gc function, stands for "gradient change". Specifies how the gradient changes, e.g. with time
#'
#' @description
#' models niche preferences by removing events (fossil occurrences) when they are outside of their preferred niche using the function `thin`
#' Combines the functions niche_def and gc ("gradient change") to determine how the taxons preferred environment changes with time. This is done by composing `niche_def` and `gc`. The result is then used as a thinning.
#' Models niche preferences by removing events (fossil occurrences) when they are outside of a preferred niche using the function `thin`.
#' Combines the functions `niche_def` and `gc` ("gradient change") to determine how the taxons' preferred environment changes with time. This is done by composing `niche_def` and `gc`. The result is then used as a thinning on the events `x`.
#'
#' @examples
#' \dontrun{
Expand All @@ -36,9 +36,13 @@ apply_niche_pref = function(x, niche_def, gc){
#' # fossil occurrences with niche preference
#' hist(foss_occ_niche, xlab = "time")
#'
#' # see also
#' vignette("event_data")
#' # for a longer example
#'
#' }
#'
#' @seealso [apply_taphonomy()] to model taphonomic effects based on the sampe principle, [thin()] for the underlying mathematical procedure
#' @seealso [apply_taphonomy()] to model taphonomic effects based on the same principle, [thin()] for the underlying mathematical procedure.

# function that returns niche preference as a function of y (typically time)
change_in_niche = function(y) niche_def(gc(y))
Expand Down
13 changes: 7 additions & 6 deletions R/apply_taphonomy.R
Original file line number Diff line number Diff line change
@@ -1,17 +1,18 @@
apply_taphonomy = function(x, pres_potential, ctc){
#' @export
#'
#' @title model taphonomy
#' @title Model taphonomic effects
#'
#' @param x numeric vector. event type data, e.g. fossil occurrences
#' @param pres_potential function. takes taphonomic conditions as input and returns the preservation potential (a number between 0 and 1)
#' @param ctc function, change in taphonomic conditions (ctc) with time.
#' @param x events, e.g. times/ages of fossil occurrences or their stratigraphic position.
#' @param pres_potential function. Takes taphonomic conditions as input and returns the preservation potential (a number between 0 and 1)
#' @param ctc function, change in taphonomic conditions (ctc) with time or stratigraphic position.
#'
#' @description
#' models taphonomy by combining the change in taphonomic conditions with the preservation potential as a function of taphonomic conditions to determine how preservation potential changes. This is then used to systematically remove (thin) the event data using `thin`
#' Models taphonomy by combining the change in taphonomic conditions with the preservation potential as a function of taphonomic conditions to determine how preservation potential changes. This is then used to systematically remove (thin) the event data using `thin`.
#'
#'
#' @seealso [apply_niche_pref()] for modeling niche preferences based on the same principle, [thin()] for the underlying mathematical procedure
#' @seealso [apply_niche_pref()] for modeling niche preferences based on the same principle, [thin()] for the underlying mathematical procedure.
#'
#'

# function that returns preservation potential as a function of input (e.g. time or position)
Expand Down
10 changes: 8 additions & 2 deletions R/cnd.R
Original file line number Diff line number Diff line change
Expand Up @@ -13,19 +13,25 @@ cnd = function(mean, sd, inc = 1, cut_neg = TRUE){
#'
#' @examples
#' \dontrun{
#' # using water depht as niche
#' # using water depth as niche
#' wd = seq(-3, 40, by = 0.5)
#' f = cnd(mean = 10, sd = 5, inc = 15, cut_neg = FALSE)
#' # 1 indicates high preference, 0 indicates low preference
#' plot(wd, f(wd), xlab = "Water depth", ylab = "Env. preference")
#' # set value at neg wd to 0 - non-terrestrial species.
#' f = cnd(mean = 10, sd = 5, inc = 15, cut_neg = TRUE)
#' plot(wd, f(wd), xlab = "Water depth", ylab = "Env. preference")
#'
#' # see also
#' vignette("event_data")
#' #for examples how to use it for niche modeling
#' }
#'
#'
#' @description
#' returns a function that defines a niche based on a capped normal distribution, i.e. a pdf of a normal distribution where all values above 1 are capped. Mathematically, this is f(x) = min( inc * pdf(x), 1)
#' Returns a function that defines a niche based on a capped normal distribution, i.e. a pdf of a normal distribution where all values above 1 are capped. Mathematically, this is f(x) = min( inc * pdf(x), 1).
#' Corresponds to the probability of collection model by Holland and Patzkowsky (2002).
#' @references * Holland, Steven M. and Patzkowsky, Mark E. 2002. "Stratigraphic variation in the timing of first and last occurrences." PALAIOS. https://doi.org/10.1669/0883-1351(2002)017<0134:SVITTO>2.0.CO;2

fa = function(x){
y = pmin(inc * stats::dnorm(x, mean, sd), rep(1, length(x)))
Expand Down
17 changes: 0 additions & 17 deletions R/cnd_niche.R

This file was deleted.

12 changes: 6 additions & 6 deletions R/ornstein_uhlenbeck.R
Original file line number Diff line number Diff line change
Expand Up @@ -3,16 +3,16 @@ ornstein_uhlenbeck = function(t, mu = 0, theta = 1, sigma = 1, y0 = 0){
#'
#' @title simulate ornstein-uhlenbeck (OU) process
#'
#' @param t times at which the process is simulated
#' @param mu numeric, long term mean
#' @param theta numeric, mean reversion speed
#' @param t times at which the process is simulated. Can be heterodistant
#' @param mu scalar, long term mean
#' @param theta scalar, mean reversion speed
#' @param sigma positive scalar, strength of randomness
#' @param y0 scalar. initial value (value of process at the first entry of t)
#' @param y0 scalar, initial value (value of process at the first entry of t)
#'
#' @description
#' simulates the Ornstein-Uhlenbeck process using the Euler-Maruyame method. The process is simulated on a scale of 0.25 * min(diff(t)) and then interpolated to the values of `t`.
#' Simulates an Ornstein-Uhlenbeck process using the Euler-Maruyama method. The process is simulated on a scale of `0.25 * min(diff(t))` and then interpolated to the values of `t`.
#'
#' @returns a list with two elements: `t` and `y`. `t` is a duplicate of the input `t`, `y` are the values of the OU process at these times. Outputs are of class `timelist` and can thus be plotted directly using `plot`, see `?plot.timelist`
#' @returns A list with two elements: `t` and `y`. `t` is a duplicate of the input `t`, `y` are the values of the OU process at these times. Output list is of class `timelist` and can thus be plotted directly using `plot`, see `?admtools::plot.timelist`
#'
#' @examples
#' \dontrun{
Expand Down
10 changes: 9 additions & 1 deletion R/p3.R
Original file line number Diff line number Diff line change
Expand Up @@ -12,19 +12,27 @@ p3 = function(rate, from, to, n = NULL){
#' @description
#' Simulates events in the interval `from` to `to` based on a Poisson point process with rate `rate`. If the parameter `n` is used, the number of fossils is conditioned to be `n`
#' In the context of paleontology, these events can be interpreted as fossil occurrences or first/last occurrences of species. In this case, the rate is the average number of fossil occurrences (resp first/last occurrences) per unit
#'
#' @returns a numeric vector with timing/location of events.
#'
#' @examples
#' \dontrun{
#' # for fossil occ.
#' x = p3(rate = 5, from = 0, to = 1) # 5 fossil occurrences per myr on avg.
#' hist(x, xlab = "Time (Myr"), ylab = "Fossil Occurrences" )
#' hist(x, xlab = "Time (Myr)", ylab = "Fossil Occurrences" )
#'
#' x = p3(rate = 3, from = 0, to = 4)
#' hist(x, main = paste0(length(x), " samples")) # no of events is random
#'
#' x = p3(rate = 3, from = 0, to = 4, n = 10)
#' hist(x, main = paste0(length(x), " samples")) # no of events is fixed to n
#'
#' # see also
#' vignette("event_data")
#' # for details on usage and applications to paleontology
#' }
#'
#' @seealso [p3_var_rate()] for the variable rate implementation

if (rate <= 0){ stop("Need strictly positive rate")}
if (from >= to) {stop("parameter \"from\" needs to be smaller than \"to\".")}
Expand Down
9 changes: 7 additions & 2 deletions R/p3_var_rate.R
Original file line number Diff line number Diff line change
Expand Up @@ -18,18 +18,23 @@ p3_var_rate = function(x, y = NULL, from = 0, to = 1, f_max = 1, n = NULL){
#' \dontrun{
#' # assuming events are fossil occurrences
#' # then rate is the avg rate of fossil occ. per unit
#' linear decrease in rate from 50 at x = 0 to 0 at x = 1
#' #linear decrease in rate from 50 at x = 0 to 0 at x = 1
#' x = c(0, 1)
#' y = c(50, 0)
#' s = p3_var_rate(x, y, f_max = 50)
#' hist(s, xlab = "Time (myr)", main = "Fossil Occurrences")
#' # conditoned to return 100 samples
#' # conditioned to return 100 samples
#' s = p3_var_rate(x, y, f_max = 50, n = 100)
#' # hand over function
#' s = p3_var_rate(x = sin, from = 0 , to = 3 * pi, n = 50)
#' hist(s) # note that negative values of f (sin) are ignored in sampling
#'
#' # see also
#' vignette("event_data")
#' # for details on usage and applications to paleontology
#' }
#'
#' @seealso [p3()] for the constant rate implementation, [rej_samp()] for the underlying random number generation.

if (from >= to){
stop("\"from\" must be smaller than \"to\".")
Expand Down
12 changes: 6 additions & 6 deletions R/random_walk.R
Original file line number Diff line number Diff line change
Expand Up @@ -3,15 +3,15 @@ random_walk = function(t, sigma = 1, mu = 0, y0 = 0){
#'
#' @title continuous time random walk
#'
#' @param t numeric vector with strictly increasing elements. Times at which the random walk is evaluated
#' @param sigma variance parameter
#' @param mu directionality parameter
#' @param y0 starting value, i.e. value of the random walk at the first entry of `t`
#' @param t numeric vector with strictly increasing elements, can be heterodistant. Times at which the random walk is evaluated
#' @param sigma positive scalar, variance parameter
#' @param mu scalar, directionality parameter
#' @param y0 scalar, starting value (value of the random walk at the first entry of `t`)
#'
#' @description
#' simulates a random walk as a Brownian drift
#' Simulates a (continuous time) random walk as a Brownian drift
#'
#' @returns a list with elements `t` and `y`. `t` is a duplicate of the input parameter and is the times at which the random walk is evaluated. `y` are the values of the random walk at said times
#' @returns A list with elements `t` and `y`. `t` is a duplicate of the input parameter and is the times at which the random walk is evaluated. `y` are the values of the random walk at said times. Output list is of class `timelist` and can thus be plotted directly using `plot`, see `?admtools::plot.timelist`
#'
#' @examples
#' \dontrun{
Expand Down
4 changes: 2 additions & 2 deletions R/rej_samp.R
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,7 @@ rej_samp = function(f, x_min, x_max, n = 1L, f_max = 1){
#' @title rejection sampling
#'
#' @description
#' rejection sampling from the (pseudo) pdf `f` in the interval between `x_min` and `x_max`. Returns `n` samples. Note that values of `f` below 0 are capped to zero
#' Rejection sampling from the (pseudo) pdf `f` in the interval between `x_min` and `x_max`. Returns `n` samples. Note that values of `f` below 0 are capped to zero
#'
#' @param f function. (pseudo) pdf from which the sample is drawn
#' @param x_min scalar. lower limit of the examined interval
Expand All @@ -18,7 +18,7 @@ rej_samp = function(f, x_min, x_max, n = 1L, f_max = 1){
#' x = rej_samp(f, 0, 3*pi, n = 100)
#' hist(x) # note that no samples are drawn where sin is negative
#' }
#'
#' @seealso [p3_var_rate()] for the derived variable rate Poisson point process implementation.
x = c()
warn = FALSE
if (f_max <= 0) {stop("`f_max` must be positive.")}
Expand Down
19 changes: 15 additions & 4 deletions R/scenarioA.R
Original file line number Diff line number Diff line change
@@ -1,4 +1,8 @@
#' @title scenario A from Hohmann et al 2024
#' @title example data, scenario A from Hohmann et al. (2024)
#'
#' @description
#' Scenario A as described in Hohmann et al. (2024), published in Hohmann et al. (2023). Contains data from a carbonate platform simulated using CarboCAT Lite (Burgess 2013, 2023)
#'
#'
#' @format A list with 6 elements:
#'
Expand All @@ -7,7 +11,14 @@
#' * `dist_from_shore` : character vector. Distance from shore in km of locations at which the observations were made
#' * `h_m` : matrix of size length(t_myr) x length(dist_from_shore). Accumulated sediment height in m at examined locations
#' * `wd_m`: matrix of size length(t_myr) x length(dist_from_shore). Water depth in m at examined locations
#' * `strat_col`: list with length(dist_from shore) elements. Represents a stratigraphic column. Each element is a list with two elements
#' * `bed_thickness_m`: numeric vector. Bed thickness in m
#' * `facies_code` : integer vector. facies code of the bed
#' * `strat_col`: list with length(dist_from shore) elements. Represents a stratigraphic column. Each element is a list with two elements:
#' * `bed_thickness_m`: numeric vector. Bed thickness in m
#' * `facies_code` : integer vector. facies code of the bed
#'
#'
#' @references
#' * Burgess, Peter. 2013. "CarboCAT: A cellular automata model of heterogeneous carbonate strata." Computers & Geosciences. <https://doi.org/10.1016/j.cageo.2011.08.026>.
#' * Burgess, Peter. 2023. "CarboCATLite v1.0.1." Zenodo. <https://doi.org/10.5281/zenodo.8402578>
#' * Hohmann, Niklas; Koelewijn, Joël R.; Burgess, Peter; Jarochowska, Emilia. 2024. "Identification of the mode of evolution in incomplete carbonate successions." BMC Ecology and Evolution, In Press. https://doi.org/10.1101/2023.12.18.572098.
#' * Hohmann, Niklas, Koelewijn, Joël R.; Burgess, Peter; Jarochowska, Emilia. 2023. “Identification of the Mode of Evolution in Incomplete Carbonate Successions - Supporting Data.” Open Science Framework. https://doi.org/10.17605/OSF.IO/ZBPWA, published under the [CC-BY 4.0](https://creativecommons.org/licenses/by/4.0/) license.
"scenarioA"
9 changes: 5 additions & 4 deletions R/stasis.R
Original file line number Diff line number Diff line change
Expand Up @@ -3,13 +3,14 @@ stasis = function(t, mean = 0, sd = 1){
#'
#' @title simulate phenotypic stasis
#'
#' @param t times at which the phenotype is determined
#' @param mean mean value
#' @param sd standard deviation
#' @param t times at which the traits are determined
#' @param mean scalar, mean trait value
#' @param sd strictly positive scalar, standard deviation of traits
#'
#' @description
#' simulates stasis as independent, normally distributed random variables with mean `mean` and standard deviation `sd`
#' Simulates stasis as independent, normally distributed random variables with mean `mean` and standard deviation `sd`
#'
#' @returns A list with two elements: `t` and `y`. `t` is a duplicate of the input `t`, `y` are the corresponding trait values. Output list is of class `timelist` and can thus be plotted directly using `plot`, see `?admtools::plot.timelist`
#' @examples
#' \dontrun{
#' library("admtools") # required for plotting of results
Expand Down
9 changes: 5 additions & 4 deletions R/thin.R
Original file line number Diff line number Diff line change
Expand Up @@ -9,17 +9,18 @@ thin = function(x, thin){
#'
#'
#' @description
#' Thins a vector of events using the function thin, meaning the probability that the ith event in x is preserved is given by _thin (x(i))_. Values of
#' `thin` below 0 and above 1 are ignored
#' Is used to model niche preferences in apply_niche_pref, where events are thinned based on how close they are to their preferred niche
#' Can also be used to model taphonomic effects. In this case, `thin` describes the preservation potential.
#' Thins a vector of events using the function thin, meaning the probability that the ith event in x is preserved is given by _thin(x(i))_. Values of
#' `thin` below 0 and above 1 are ignored.
#' Is used to model niche preferences in `apply_niche_pref` and taphonomic effects in `apply_taphonomy`.
#'
#' @examples
#' \dontrun{
#' x = p3(rate = 100, from = 0, to = 3 * pi) # simulate Poisson point process
#' y = thin(x, sin)
#' hist(y) # not how negative values of sin are treated as 0
#' yy = thin(x, function(x) 5 * sin(x))
#' hist(yy) # note how values of 5 * sin above 1 are not affecting the thinning
#' }
#'
#' @seealso [apply_niche_pref()] and [apply_taphonomy()] for use cases with biological meaning

Expand Down
24 changes: 24 additions & 0 deletions man/StratPal-package.Rd

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

Loading
Loading