Merge pull request #6 from MindTheGap-ERC/dev

Dev

.Rbuildignore

@@ -7,3 +7,5 @@
 ^pkgdown$
 
 ^CONTRIBUTING.md$
+^doc$
+^Meta$

.gitignore

@@ -1,3 +1,5 @@
 .Rproj.user
 inst/doc
 docs
+/doc/
+/Meta/

NAMESPACE

@@ -1,6 +1,7 @@
 # Generated by roxygen2: do not edit by hand
 
 export(apply_niche_pref)
+export(apply_taphonomy)
 export(cnd)
 export(ornstein_uhlenbeck)
 export(p3)

R/apply_niche_pref.R

@@ -38,7 +38,11 @@ apply_niche_pref = function(x, niche_def, gc){
     #'
     #' }
     #'
-  niche_with_time = function(t) niche_def(gc(t))
-  r = thin(x, niche_with_time)
+  #' @seealso [apply_taphonomy()] to model taphonomic effects based on the sampe 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))
+  # thin events based on niche preference
+  r = thin(x, change_in_niche)
   return(r)
   }

R/apply_taphonomy.R

@@ -0,0 +1,22 @@
+apply_taphonomy = function(x, pres_potential, ctc){
+  #' @export
+  #'
+  #' @title model taphonomy
+  #'
+  #' @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.
+  #'
+  #' @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`
+    #'
+    #'
+  #' @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)
+  change_pres_pot = function(y) pres_potential(ctc(y))
+  # thin events based on preservation potential
+  r = thin(x, change_pres_pot)
+  return(r)
+}

R/ornstein_uhlenbeck.R

@@ -10,7 +10,7 @@ ornstein_uhlenbeck = function(t, mu = 0, theta = 1, sigma = 1, y0 = 0){
   #' @param y0 scalar. initial value (value of process at the first entry of t)
   #'
   #' @description
-    #' simulates ornstein-uhlenbeck 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 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`.
     #'
   #' @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`
   #'

R/stasis.R

@@ -8,7 +8,7 @@ stasis = function(t, mean = 0, sd = 1){
   #' @param sd standard deviation
   #'
   #' @description
-    #' simulates stasis as independent, normally distributed random variables with mean `mean` and standard deviatin `sd`
+    #' simulates stasis as independent, normally distributed random variables with mean `mean` and standard deviation `sd`
     #'
   #' @examples
     #' \dontrun{

R/thin.R

@@ -11,6 +11,8 @@ 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.
     #'
   #' @examples
     #' x = p3(rate = 100, from = 0, to = 3 * pi) # simulate Poisson point process
@@ -19,10 +21,11 @@ thin = function(x, thin){
     #' 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
 
   p = pmax(pmin(thin(x), rep(1, length(x))), rep(0, length(x))) # cut off at 0 and 1
-  ind = stats::rbinom(n = length(x), size = 1, prob = p)
-  y = x[as.logical(ind)]
+  ind = stats::rbinom(n = length(x), size = 1, prob = p) # determine if preserved or not
+  y = x[as.logical(ind)] # select preserved events
   return(y)
 
 }

vignettes/advanced_functionality.Rmd

@@ -18,4 +18,65 @@ knitr::opts_chunk$set(
 library(StratPal)
 ```
 
-TBA
+Here we go through how the presented workflows can be generalized and expanded upon.
+
+## Different forward models
+
+In the examples, we used age-depth models from a carbonate platform simulated with CarboCAT lite. This can easily be expanded to any other sedimentary forward model (be it marine or terrestrial, siliciclastic, mixed, or carbonate system).
+
+For this, age-depth information needs to be extracted from the forward model. The vectors of elapsed model time vs. accumulated height can then be handed over to `tp_to_adm` to define age-depth model objects. These can then be used in the pipelines as in the examples.
+
+If the forward model included erosion, you need to first account for time intervals where sediment is first deposited and later removed. For this, first define a sediment accumulation curve using `tp_to_sac`, and pass it to `sac_to_adm` to turn it into an age-depth model. You can then use the resulting age-depth model for your analysis pipeline.
+
+## Empirical age-depth models
+
+Of course you can use empirical age-depth models in your pipeline. You can create them from tie points using `tp_to_adm` just as in the case with forward models.
+
+The workflow in `StratPal` does not work with `multiadm` objects as defined by `admtools`. These objects store multiple age-depth models to account for uncertainty. You can however reduce their complexity, e.g. via `admtools::median_adm`, and use the resulting age-depth model for the modeling pipeline.
+
+## Other ecological gradients and niche definitions
+
+In the examples we used water depth as ecological gradient. This was because water depth was stored by the forward model used. The approach for niche modeling used in `StratPal` works for all types of niches. Potentially relevant niches where information can be extracted from forward models are e.g. substrate consistency, temperature, water energy, and many others.
+
+We used a simple niche definition based on a probability density function of a normal distribution. This can be expanded to arbitrary niche definitions along a gradient.
+
+For niche modeling, two components are needed:
+
+1.  A function `niche_def` . It takes as input values of the gradient, and returns numbers returns numbers between 0 and 1 as output (where 1 corresponds to the taxon being in its preferred niche, and 0 to it being completely out of its preferred niche). It must take vectorized input, and return a vector of equal length
+2.  A function `gc` that describes how the gradient changes with time. It must take vectorized inputs, and return a vector of equal length.
+
+As long as these conditions are met, niche preferences along any gradient can be modeled using `apply_niche_pref`.
+
+## Different models of phenotypic evolution
+
+Other types of phenotypic evolution in the time domain can easily be incorporated. They should stick to the following conventions to be seamlessly integrated:
+
+-   They should take a vector of times as first argument
+
+-   The implementation should be for a continuous time version of the mode of interest, as the time domain is generally sampled at irregular times (see Hohmann et al. 2024, methods section, for details)
+
+-   They should return a list with two elements: One named `t` that is a duplicate of the first argument handed to the function. The second should be named `y` and should contain the simulated trait values.
+
+-   The output list `l` should be assigned both the class "list" and "timelist" via the command `class(l) = c("timelist", "list")` , see `random_walk` as example
+
+The last two points make sure plotting via `plot.timelist` and `plot.stratlist` works seamlessly.
+
+## Different models of event type data
+
+We used (constant and variable rate) Poisson point processes to simulate event type data. Any model of event-type data can be used (e.g. with temporal correlation), as long as it returns a vector with the timing of the events
+
+## Taphonomy for event type data
+
+Taphonomic effects can be incorporated using using the `apply_taphonomy` function. This is based on the same principle as niche modeling, and makes use of the `thin` function, but was not shown in the examples.
+
+## Transform different types of data
+
+`strat_to_time` and `time_to_strat` are generic functions of the `admtools` package that can transform different types of data. Currently they transform phylogenetic trees (`phylo` objects), lists (`list` class), and numeric data such as the event type data (class `numeric`).
+
+In general, any type temporal (resp. stratigraphic) data can be transformed. For this, you need define an `S3` class for the object you would like to transform, and then implement the transformation of this object for `time_to_strat` and `strat_to_time`.
+
+If you would like to add transformations for new S3 classes to the `admtools` package, please see the contribution guidelines of the `admtools` package for details on how to contribute code.
+
+## References
+
+\* Niklas Hohmann, Joël R Koelewijn, Peter Burgess, Emilia Jarochowska. Identification of the mode of evolution in incomplete carbonate successions.bioRxiv 2023.12.18.572098; doi: [10.1101/2023.12.18.572098]([https://doi.org/10.1101/2023.12.18.572098),](https://doi.org/10.1101/2023.12.18.572098),) Data published under a [CC-BY 4.0]([https://creativecommons.org/licenses/by/4.0/)](https://creativecommons.org/licenses/by/4.0/)) license.