From b96e18d5e0af253f259725edfa191eeb834c68ca Mon Sep 17 00:00:00 2001
From: Niklas Hohmann <67792281+NiklasHohmann@users.noreply.github.com>
Date: Fri, 19 Jul 2024 19:11:33 +0200
Subject: [PATCH 1/5] typos
---
vignettes/phenotypic_evolution.Rmd | 31 +++++++++++++++---------------
1 file changed, 16 insertions(+), 15 deletions(-)
diff --git a/vignettes/phenotypic_evolution.Rmd b/vignettes/phenotypic_evolution.Rmd
index 54f372b..e78d21c 100644
--- a/vignettes/phenotypic_evolution.Rmd
+++ b/vignettes/phenotypic_evolution.Rmd
@@ -42,7 +42,7 @@ Use `vignette("advanced_functionality")` for guidelines to implement your own mo
You can visualize the different modes of evolution using the following pipeline:
```{r}
-seq(0, 1, by = 0.01) |> # times of simulation in myr. simulate over 1 Myr years with 10 kyr resolution
+seq(0, 1, by = 0.01) |> # times of simulation in Myr. Simulate over 1 Myr years with 10 kyr resolution
random_walk(sigma = 1, mu = 3) |> # simulate random walk with increasing trait values
plot(type = "l", # plot results
xlab = "Time [Myr]",
@@ -53,15 +53,15 @@ seq(0, 1, by = 0.01) |> # times of simulation in myr. simulate over
## Stratigraphic domain
-We are interested in how phenotypic evolution is preserved in the stratigraphic record, i.e. how is trait evolution within a lineage preserved at a specific location in the stratigraphic record. Here we develop the modeling pipelines to answer this question.
+We are interested in how phenotypic evolution is preserved in the stratigraphic record, i.e. how trait evolution within a lineage is preserved at a specific location in the stratigraphic record. Here we develop the modeling pipelines to answer this question.
### Age-depth models
-As an example, we compare the the preservation of trait evolution 2 km and 12 km offshore in the carbonate platform in scenario A. See `?scenarioA` and `vignette("StratPal")` for details on scenario A.
+As an example, we compare the preservation of trait evolution 2 km and 12 km offshore in the carbonate platform in scenario A. See `?scenarioA` and `vignette("StratPal")` for details on scenario A.
First define the age-depth models:
-```{r, figures-side, fig.show="hold", out.width="50%", fig.align='left'}
+```{r, figures-side, fig.show="hold", out.width="50%"}
adm_2km = tp_to_adm(t = scenarioA$t_myr, # 2 km from shore
h = scenarioA$h_m[,"2km"],
T_unit = "Myr",
@@ -72,7 +72,6 @@ adm_12km = tp_to_adm(t = scenarioA$t_myr, # 12 km from shore
T_unit = "Myr",
L_unit = "m")
-
plot(adm_2km, # plot age-depth model 2 km from shore
lwd_acc = 2, # plot thicker lines for intervals with sediment accumulation (lwd = line width)
lty_destr = 0) # don't plot destructive intervals/gaps (lty = line type)
@@ -100,8 +99,8 @@ seq(from = min_time(adm_2km), to = max_time(adm_2km), by = 0.01) |> # sample eve
plot(type = "l", # plot
orientation = "lr",
xlab = paste0("Stratigraphic height [", get_L_unit(adm_2km), "]"),
- ylab = "Trait value")
-
+ ylab = "Trait value",
+ main = "Trait evolution 2 km from shore")
```
This is what a biased random walk in the time domain would look like if it was observed 2 km offshore in the simulated carbonate platform. The large jumps in traits correspond to long hiatuses caused by prolonged drops in relative sea level (Hohmann et al. 2024). Compare this figure to the random walk in the time domain (without stratigraphic distortions).
@@ -115,10 +114,11 @@ seq(from = min_time(adm_12km), to = max_time(adm_12km), by = 0.01) |> # sample e
plot(type = "l", # plot results
orientation = "lr",
xlab = paste0("Stratigraphic height [", get_L_unit(adm_12km), "]"),
- ylab = "Trait value")
+ ylab = "Trait value",
+ main = "Trait evolution 12 km from shore")
```
-You can also see two jumps, but the first is at a different location compared to the section 2 km offshore. Looking at the age-depth models, we can see why this is: 2 km offshore, the jumps are caused by prolonged hiatuses that are associated with the massive drops in sea level around 0.5 Myr and 1.5 Myr. 12 km offshore, the age-depth model shows fewer gaps, but prolonged intervals with low sedimentation rates at the beginning of the simulation and at around 1.5 Myr. This leads to intervals where the rate phenotypic evolution appears to be accelerated because of stratigraphic condensation (Hohmann 2021). You can see that not all artefactual jumps in traits observable in the stratigraphic record are caused by gaps.
+You can also see two jumps, but the first is at a different location compared to the section 2 km offshore. Looking at the age-depth models, we can see why this is: 2 km offshore, the jumps are caused by prolonged hiatuses that are associated with the massive drops in sea level around 0.5 Myr and 1.5 Myr. 12 km offshore, the age-depth model shows fewer gaps, but prolonged intervals with low sedimentation rates at the beginning of the simulation and at around 1.5 Myr. This leads to intervals where the rate phenotypic evolution appears to be accelerated because of (sedimentological) condensation (Hohmann 2021). You can see that not all artefactual jumps in traits observable in the stratigraphic record are caused by gaps.
**Task:** How does this effect vary between different modes of evolution (with different parameter choices), and at different locations in the platform (e.g., along an onshore-offshore gradient)? Can you make any general statements about where you see the strongest effects?
@@ -130,11 +130,11 @@ However, there is a small imperfection. Because we simulate the lineage every 10
To fix this, we want to prescribe where in the stratigraphic domain samples are taken. This allows us to examine the effect that different sampling strategies have on our perception of trait evolution. We do this in the following steps:
-1. Determine sampling locations in the stratigraphic domain
-2. Calculate what times correspond to said sampling locations via `strat_to_time`
-3. Simulate a lineage at said times
-4. Transform the simulated trait values back into the stratigraphic domain using `time_to_strat`
-5. Plot the result
+1. Determine sampling locations in the stratigraphic domain.
+2. Calculate what times correspond to said sampling locations via `strat_to_time`.
+3. Simulate a lineage at said times.
+4. Transform the simulated trait values back into the stratigraphic domain using `time_to_strat`.
+5. Plot the result.
Let's assume we take a sample every 2 meters. Then our five steps above result in the following modeling pipeline:
@@ -151,7 +151,8 @@ sampling_loc_m |> # sampling locations
plot(orientation = "lr", # plot stratigrahic data
type = "l",
ylab = "Trait Value",
- xlab = paste0("Stratigraphic height [", get_L_unit(adm_2km), "]"))
+ xlab = paste0("Stratigraphic height [", get_L_unit(adm_2km), "]"),
+ main = "Trait evolution 2 km from shore")
```
There we have it! This is what the lineage would look like if it was sampled every 2 meters in a sections 2 km from shore.
From 5179e47e2831c63910d4febb4257e96a220a69da Mon Sep 17 00:00:00 2001
From: Niklas Hohmann <67792281+NiklasHohmann@users.noreply.github.com>
Date: Sat, 20 Jul 2024 08:44:22 +0200
Subject: [PATCH 2/5] add links to refs
---
vignettes/phenotypic_evolution.Rmd | 12 ++++++------
1 file changed, 6 insertions(+), 6 deletions(-)
diff --git a/vignettes/phenotypic_evolution.Rmd b/vignettes/phenotypic_evolution.Rmd
index e78d21c..27081d4 100644
--- a/vignettes/phenotypic_evolution.Rmd
+++ b/vignettes/phenotypic_evolution.Rmd
@@ -29,11 +29,11 @@ library(admtools)
The `StratPal` package provides three functions for simulating different modes of evolution in the time domain:
-- `stasis` simulates evolutionary stasis as independent, normally distributed random variables with mean `mean` and standard deviation `sd` (Hunt 2006).
+- `stasis` simulates evolutionary stasis as independent, normally distributed random variables with mean `mean` and standard deviation `sd` ([Hunt 2006](#References)).
-- `random_walk` simulates trait evolution following a (potentially biased) random walk with variability `sigma`, directionality `mu`, and initial trait value `y0` (e.g. Bookstein 1987). Setting `mu` to a negative (positive) value will make the value of the random walk decrease (increase) over time, `sigma` determines the effect of randomness on the trait values.
+- `random_walk` simulates trait evolution following a (potentially biased) random walk with variability `sigma`, directionality `mu`, and initial trait value `y0` (e.g. [Bookstein 1987](#References)). Setting `mu` to a negative (positive) value will make the value of the random walk decrease (increase) over time, `sigma` determines the effect of randomness on the trait values.
-- `ornstein_uhlenbeck` corresponds to convergence to a phenotypic optimum, where `mu` is the optimal value (long term mean), `theta` determines how fast `mu` is approached, `sigma` the noise level, and `y0` the initial trait value (Lande 1976).
+- `ornstein_uhlenbeck` corresponds to convergence to a phenotypic optimum, where `mu` is the optimal value (long term mean), `theta` determines how fast `mu` is approached, `sigma` the noise level, and `y0` the initial trait value ([Lande 1976](#References)).
Use `vignette("advanced_functionality")` for guidelines to implement your own modes of evolution.
@@ -103,7 +103,7 @@ seq(from = min_time(adm_2km), to = max_time(adm_2km), by = 0.01) |> # sample eve
main = "Trait evolution 2 km from shore")
```
-This is what a biased random walk in the time domain would look like if it was observed 2 km offshore in the simulated carbonate platform. The large jumps in traits correspond to long hiatuses caused by prolonged drops in relative sea level (Hohmann et al. 2024). Compare this figure to the random walk in the time domain (without stratigraphic distortions).
+This is what a biased random walk in the time domain would look like if it was observed 2 km offshore in the simulated carbonate platform. The large jumps in traits correspond to long hiatuses caused by prolonged drops in relative sea level ([Hohmann et al. 2024](#References)). Compare this figure to the random walk in the time domain (without stratigraphic distortions).
12 km offshore, the preservation is very different:
@@ -118,7 +118,7 @@ seq(from = min_time(adm_12km), to = max_time(adm_12km), by = 0.01) |> # sample e
main = "Trait evolution 12 km from shore")
```
-You can also see two jumps, but the first is at a different location compared to the section 2 km offshore. Looking at the age-depth models, we can see why this is: 2 km offshore, the jumps are caused by prolonged hiatuses that are associated with the massive drops in sea level around 0.5 Myr and 1.5 Myr. 12 km offshore, the age-depth model shows fewer gaps, but prolonged intervals with low sedimentation rates at the beginning of the simulation and at around 1.5 Myr. This leads to intervals where the rate phenotypic evolution appears to be accelerated because of (sedimentological) condensation (Hohmann 2021). You can see that not all artefactual jumps in traits observable in the stratigraphic record are caused by gaps.
+You can also see two jumps, but the first is at a different location compared to the section 2 km offshore. Looking at the age-depth models, we can see why this is: 2 km offshore, the jumps are caused by prolonged hiatuses that are associated with the massive drops in sea level around 0.5 Myr and 1.5 Myr. 12 km offshore, the age-depth model shows fewer gaps, but prolonged intervals with low sedimentation rates at the beginning of the simulation and at around 1.5 Myr. This leads to intervals where the rate phenotypic evolution appears to be accelerated because of (sedimentological) condensation ([Hohmann 2021](#References)). You can see that not all artefactual jumps in traits observable in the stratigraphic record are caused by gaps.
**Task:** How does this effect vary between different modes of evolution (with different parameter choices), and at different locations in the platform (e.g., along an onshore-offshore gradient)? Can you make any general statements about where you see the strongest effects?
@@ -177,7 +177,7 @@ vignette("advanced_functionality")
for details on how to expand on the modeling pipelines described here, or explore the vignette online under [mindthegap-erc.github.io/StratPal/articles/advanced_functionality](https://mindthegap-erc.github.io/StratPal/articles/advanced_functionality.html).
-## References
+## References {#References}
* Bookstein, Fred L. 1987. "Random walk and the existence of evolutionary rates." Paleobiology. https://doi.org/10.1017/S0094837300009039.
From 774016d9cfd9d6dc015319f91b864797262d4fba Mon Sep 17 00:00:00 2001
From: Niklas Hohmann <67792281+NiklasHohmann@users.noreply.github.com>
Date: Sat, 20 Jul 2024 09:45:55 +0200
Subject: [PATCH 3/5] fix ref
---
vignettes/phenotypic_evolution.Rmd | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
diff --git a/vignettes/phenotypic_evolution.Rmd b/vignettes/phenotypic_evolution.Rmd
index 27081d4..095c2ef 100644
--- a/vignettes/phenotypic_evolution.Rmd
+++ b/vignettes/phenotypic_evolution.Rmd
@@ -57,7 +57,7 @@ We are interested in how phenotypic evolution is preserved in the stratigraphic
### Age-depth models
-As an example, we compare the preservation of trait evolution 2 km and 12 km offshore in the carbonate platform in scenario A. See `?scenarioA` and `vignette("StratPal")` for details on scenario A.
+As an example, we compare the preservation of trait evolution 2 km and 12 km offshore in the carbonate platform in scenario A from [Hohmann et al. 2024](#References). See `?scenarioA` and `vignette("StratPal")` for details on scenario A.
First define the age-depth models:
From e31232efaa014b3bec31bafea2a4a7c37ed78a83 Mon Sep 17 00:00:00 2001
From: Niklas Hohmann <67792281+NiklasHohmann@users.noreply.github.com>
Date: Sat, 20 Jul 2024 17:24:10 +0200
Subject: [PATCH 4/5] expand vignettes
---
vignettes/StratPal.Rmd | 133 ++++++++------
vignettes/advanced_functionality.Rmd | 53 +++---
vignettes/event_data.Rmd | 254 ++++++++++++++++-----------
vignettes/phenotypic_evolution.Rmd | 2 +-
4 files changed, 254 insertions(+), 188 deletions(-)
diff --git a/vignettes/StratPal.Rmd b/vignettes/StratPal.Rmd
index 6c1a2d0..af25480 100644
--- a/vignettes/StratPal.Rmd
+++ b/vignettes/StratPal.Rmd
@@ -14,42 +14,30 @@ knitr::opts_chunk$set(
)
```
-```{r setup}
+```{r setup, include=FALSE}
library(StratPal)
library(admtools)
```
## Overview
-Welcome to the `StratPal` package. This vignette provides an overview of the structure of the package and preliminaries needed to efficiently use it. We go through _installation_, _dependencies_, provide an overview of the available _example data_, _piping_, and working with _age-depth models_.
+Welcome to the `StratPal` package. This vignette provides an overview of the structure of the package and preliminaries needed to efficiently use it. We go through [*installation*](#installation), [*dependencies*](#dependencies), provide an overview of the available [*example data*](#example-data), [*piping*](#piping), and working with [*age-depth models*](#adm).
-Throughout the vignettes there are several **tasks** that you can solve. They are not mandatory for understanding the functionality of the package. Their aim is to explore the available models and data and develop an intuition for stratigraphic paleobiology.
+If you want to skip the introduction, you can also [get started](#getting-started) right away.
-If you want to skip the introduction, go to
+Throughout the vignettes there are several **tasks** that you can solve. They are not required for understanding the functionality of the package. Their aim is to explore the available models and data and develop an intuition for stratigraphic paleobiology.
-```{r, eval=FALSE}
-vignette("phenotypic_evolution")
-```
-
-for details on how to model stratigraphic paleobiology of phenotypic evolution, or explore the vignette online on the package webpage here: [mindthegap-erc.github.io/StratPal/articles/phenotypic_evolution](https://mindthegap-erc.github.io/StratPal/articles/phenotypic_evolution.html). Go to
-
-```{r, eval=FALSE}
-vignette("event_data")
-```
-
-for details on how to model stratigraphic paleobiology of event data such as individual fossils and first/last occurrences of taxa, or explore the vignette online on the package website here: [mindthegap-erc.github.io/StratPal/articles/event_data](https://mindthegap-erc.github.io/StratPal/articles/event_data.html).
-
-## Installation
+## Installation {#installation}
To install the `StratPal` package, first install the `remotes` package by running
-```R
+``` r
install.packages("remotes")
```
in the R console. Then, run
-```R
+``` r
remotes::install_github(repo = "MindTheGap-ERC/StratPal",
build_vignettes = TRUE,
ref = "HEAD",
@@ -58,32 +46,48 @@ remotes::install_github(repo = "MindTheGap-ERC/StratPal",
to install the latest stable version of the package and all its dependencies.
-## Dependencies
+## Dependencies {#dependencies}
-The `StratPal` heavily relies on the age-depth modelling tools provided by the `admtools` package, which is _automatically_ installed when you install `StratPal`.
+The `StratPal` heavily relies on the age-depth modelling tools provided by the `admtools` package, which is *automatically* installed when you install `StratPal`.
To use its functionality, you need to run
-```{r}
+```{r, eval=FALSE}
library(admtools)
```
-
-before running any of the provided examples. Below we provide a brief overview of the functionality of the `admtools` package we need. If you want more information, you can browse through the package vignettes using
+before running any of the examples. Below we provide a brief overview of the functionality of the `admtools` package we need. If you want more information, you can browse through the package vignettes using
```{r, eval=FALSE}
browseVignettes(package = "admtools") # opens in Browser
```
-or by visiting the package website at https://mindthegap-erc.github.io/admtools/.
+or by visiting the package website at .
-## Example data
+## Example data {#example-data}
-`StratPal` comes with some example data. This data is taken from Hohmann et al. (2024), who simulated a carbonate platform over 2 Myr. More details can be found via `?scenarionA`.
+`StratPal` comes with synthetic example data, which contains model parameters and outputs of a carbonate platform simulated using CarboCAT Lite ([Burgess (2013)](#References), [Burgess (2023)](#References)).
+The data is stored in the variable `scenarioA`, which is automatically available when the package is loaded. The structure of the data is described under `?scenarioA`.
-## Piping
+The data is taken from scenario A in [Hohmann et al. (2023)](#References), the modeling procedure is described in detail in [Hohmann et al. 2024](#References), see therein for a chronostratigraphic diagram and a transect through the carbonate platform.
-In the `StratPal` package and its vignettes, we make heavy use of the base R _pipe_ operator `|>`. While this is not required to run the package, it simplifies the code and makes the underlying logic of a modeling _pipeline_ clearer.
+The data contains information on the eustatic sea level curve used for the model run, elapsed model time, as well information on accumulated sediment thickness, water depth, and bed thicknesses and facies at locations 2, 4, 6, 8, 10, and 12 km from shore in the simulated carbonate platform.
+
+As an example, here is the eustatic sea level curve used for the model run:
+
+```{r}
+plot(x = scenarioA$t_myr,
+ y = scenarioA$sl_m,
+ type = "l",
+ xlab = "Time [Myr]",
+ ylab = "Eustatic sea level [m]",
+ main = "Sea level curve used as model input")
+```
+
+
+## Piping {#piping}
+
+In the `StratPal` package and its vignettes, we make heavy use of the base R *pipe* operator `|>`. While this is not required to run the package, it simplifies the code and makes the underlying logic of a modeling *pipeline* clearer.
This functionality is available from R version 4.2 on.
@@ -97,7 +101,7 @@ l = random_walk(t) # simulate the random walk
plot(l, type = "l") # line plot of the results
```
-This code does the job, but it has some flaws: We introduced a lot of intermediate variables, which makes it hard to trace the logic of what we are doing: plotting a random walk. Using the pipe operator `|>` we can clarify the logic and simplify the code:
+This code does the job, but it has some flaws: We introduced a lot of intermediate variables, which makes it hard to trace the logic of what we are trying to achieve: plotting a random walk. Using the pipe operator `|>` we can clarify the logic and simplify the code:
```{r}
seq(0, 1, by = 0.01) |> # define times of simulation
@@ -106,7 +110,7 @@ seq(0, 1, by = 0.01) |> # define times of simulation
```
You see that the code does the same thing: it plots a random walk, but it does so in one step without intermediate variables by chaining together the commands using the pipe `|>`. This becomes a powerful tool once we combine more and more components into longer modeling pipelines. It also makes the code readable, as you can simply read it from left to right, without having to track any intermediate steps.
-
+
Semantically, you can read the `|>` as "take the data on the left of `|>` and use it as the first argument in the function to the right of `|>`.
### Advanced usage
@@ -119,7 +123,7 @@ seq(0, 1, by = 0.1) |>
quantile(x = runif(100), p = _) # pass left hand side to the p argument
```
-## Age-depth models
+## Age-depth models {#adm}
The `StratPal` package makes heavy use of age-depth modeling functionality from the `admtools` package. Here we go through some basics of dealing with age-depth models. For more details on available functionality you can browse through the package vignettes using
@@ -127,7 +131,7 @@ The `StratPal` package makes heavy use of age-depth modeling functionality from
browseVignettes(package = "admtools") # opens in Browser
```
-or visit the package website at https://mindthegap-erc.github.io/admtools/.
+or visit the package website at .
To get started, first load the package using
@@ -135,36 +139,38 @@ To get started, first load the package using
library("admtools")
```
-The `StratPal` package comes with some example data fro age-depth models, which are available in the variable `scenarioA`, see `?scenarioA` for details on it content.
+The `StratPal` package comes with some example data for age-depth models stored in the `scenarioA` variable, see the section on [example data](#example-data) or `?scenarioA` for details.
### Defining age-depth models
Let's start with defining an age-depth model. This can be done with `tp_to_adm` (tie points to age-depth model):
```{r}
-t = scenarioA$t_myr # extract time tie points
-h = scenarioA$h_m[,"2km"] # get height tie points 2 km offshore in the carbonate plarform of scenario A from Hohmann et al. 2024
+t = scenarioA$t_myr # extract time tie points
+h = scenarioA$h_m[,"2km"] # get height tie points 2 km offshore in scenario A
# define age-depth model
# h[i] is the stratigraphic position at time t[i]
-adm = tp_to_adm(t = t, # tie points in time
- h = h, # tie points at height
- T_unit = "Myr", # time unit of adm
- L_unit = "m") # Length unit
+adm = tp_to_adm(t = t, # tie points in time
+ h = h, # tie points at height
+ T_unit = "Myr", # add time unit
+ L_unit = "m") # add length unit
```
+The nice thing about constructing age-depth models this way is that there is plenty of functionality available to extract data from age-depth models, plot them, or transform data using them.
### Plotting
-Now you can plot the adm using the basic `plot` command:
+Now you can plot the age-depth model using the basic `plot` command:
```{r}
# plot age-depth model, see ?plot.adm for details
plot(adm,
- lwd_acc = 2, # plot thicker lines for accumulative intervals (lwd = line width)
- lty_destr = 0) # don't plot destructive intervals (lty = line type)
-T_axis_lab() # add time axis label
-L_axis_lab() # add length axis label
+ lwd_acc = 2, # plot thicker lines for intervals with sediment accumulation (lwd = line width)
+ lty_destr = 0) # don't plot destructive intervals/gaps (lty = line type)
+T_axis_lab() # add time axis label
+L_axis_lab() # add length axis label
+title("Age-depth model 2 km from shore")
```
### Extracting information
@@ -177,6 +183,7 @@ get_total_thickness(adm) # sediment accumulated
get_completeness(adm) # stratigraphic completeness
summary(adm) # some summary statistics
```
+
We can now use the pipe operator to do some first analysis of the age-depth model
```{r}
@@ -187,30 +194,46 @@ adm |>
hist( xlab = paste("Hiatus duration", "[", get_T_unit(adm),"]"))
```
-Here it becomes clear that there are a 8 shorter hiatuses (below 100 kyr) and 2 long hiatuses with a duration of more than 500 kyr.
+You can see that there are a 8 shorter hiatuses (below 100 kyr) and 2 long hiatuses with a duration of more than 500 kyr. To get a detailed list with information on hiatuses you can use `get_hiat_list`.
### Transforming data
-Given a stratigraphic position, an age-depth model can tell us how old that positions is. Conversely, if we know the timing of an event, an age-depth model can tell us at what stratigraphic position that event will occur. This can be used to transfrom all types of data from the time domain to the stratigraphic domain and vice versa.
+Given a stratigraphic position, an age-depth model can tell us how old that positions is. Conversely, if we know the timing of an event, an age-depth model can tell us at what stratigraphic position that event will occur. This can be used to transform all types of data from the time domain to the stratigraphic domain and vice versa.
+
+In `admtools`, the transformation of data is done by the functions `time_to_strat` (for transforming temporal data into stratigraphic data) and `strat_to_time` (for transforming stratigrahic data into temporal data). Details on how this is done and what types of data can be transformed can be found in the [vignette of the `admtools` package](https://mindthegap-erc.github.io/admtools/articles/admtools.html), for applications to stratigraphic paleobiology see the vignettes linked below under ["Getting started"](#getting-started).
+
+## Getting started {#getting-started}
+
+With the preliminaries out of the way, you can go to
-In `admtools`, the transformation of data is done by the functions `time_to_strat` (for transforming temporal data into stratigraphic data) and `strat_to_time` (for transforming stratigrahic data into temporal data). Details on how this is done and what types of data can be transformed can be found in the other vignettes (see section below).
+```{r, eval=FALSE}
+vignette("phenotypic_evolution")
+```
-## Getting started
+for details on how to model stratigraphic paleobiology of phenotypic evolution, or explore the vignette online under [mindthegap-erc.github.io/StratPal/articles/phenotypic_evolution](https://mindthegap-erc.github.io/StratPal/articles/phenotypic_evolution.html).
-With the preliminaries out of the way, you can go to go to
+Go to
```{r, eval=FALSE}
vignette("event_data")
```
-for details on how to model stratigraphic paleobiology of events such as fossil occurrences and first/last occurrences. Go to
+for details on how to model stratigraphic paleobiology of event data such as individual fossils and first/last occurrences of taxa, or explore the vignette online under [mindthegap-erc.github.io/StratPal/articles/event_data](https://mindthegap-erc.github.io/StratPal/articles/event_data.html).
+
+See also
```{r, eval=FALSE}
-vignette("phenotypic_evolution")
+vignette("advanced_functionality")
```
-for details on how to model stratigraphic paleobiology of trait evolution.
+for details on how to expand on the modeling pipelines described here, or explore the vignette online under [mindthegap-erc.github.io/StratPal/articles/advanced_functionality](https://mindthegap-erc.github.io/StratPal/articles/advanced_functionality.html).
+
+## References {#References}
+
+- Burgess, Peter. 2013. "CarboCAT: A cellular automata model of heterogeneous carbonate strata." Computers & Geosciences. .
+
+- Burgess, Peter. 2023. "CarboCATLite v1.0.1." Zenodo.
-## References
+- 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. .
-* 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), Data published under a [CC-BY 4.0](https://creativecommons.org/licenses/by/4.0/) license.
+- 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. , published under the [CC-BY 4.0](https://creativecommons.org/licenses/by/4.0/) license.
diff --git a/vignettes/advanced_functionality.Rmd b/vignettes/advanced_functionality.Rmd
index fc50b3e..d7b01ad 100644
--- a/vignettes/advanced_functionality.Rmd
+++ b/vignettes/advanced_functionality.Rmd
@@ -14,69 +14,72 @@ knitr::opts_chunk$set(
)
```
-```{r setup}
-library(StratPal)
-```
-
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).
+In the examples, we used age-depth models from a carbonate platform simulated with CarboCAT lite ([Burgess 2013](#References), [Burgess 2023](#References)). 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.
+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 `admtools::tp_to_sac`, and pass it to `admtools::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.
+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.
## 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.
+In the examples we used water depth as ecological gradient. The approach for niche modeling used in `StratPal` works for all types of niches and gradient. Potentially relevant gradients that can be extracted from forward models are for example substrate consistency, temperature, or water energy.
+
+We used a simple niche definition based on a probability density function of a normal distribution implemented in `cnd`. This can be expanded to arbitrary niche definitions along a gradient. For integration with `apply_niche_pref`, niche definitions must meet the following criteria:
-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.
+1. They are functions that take gradient values as input and return numbers between 0 and 1 as outputs.
+2. They are vectorized, meaning when handed a vector, they return a vector of equal length.
-For niche modeling, two components are needed:
+Similar criteria hold for the definitions of gradient change `gc`:
-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.
+1. `gc` must be a function that takes time (or position) as input, and returns gradient values as outputs.
+2. The function must be vectorized, meaning when handed a vector, it must 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`.
+As long as these conditions for `gc` and `niche_def` are met, arbitrary 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:
+Other types of phenotypic evolution in the time domain can easily be incorporated. For seamless integration, they need to meet the following criteria:
-- They should take a vector of times as first argument
+- They 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)
+- The implementation is for a continuous time version of the mode of interest, as the time domain is generally sampled at irregular times (see methods section in [Hohmann et al. (2024)](#References) for details). More specifically, the implementation must be able to deal with heterodistant sampling in time.
-- 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.
+- They return a list with two elements: One named `t` that is a duplicate of the first argument handed to the function. The second one should be named `y` and 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 output list `l` should be assigned both the class "list" and "timelist" via the command `class(l) = c("timelist", "list")` , see the source code of `random_walk` for an example
-The last two points make sure plotting via `plot.timelist` and `plot.stratlist` works seamlessly.
+The last two points make sure plotting via `admtools::plot.timelist` and `admtools::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
+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., one with temporal correlation), as long as it returns a vector with the timing/position 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.
+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. For `pres_potential` (resp. `ctc`), the same logical constraints apply as to `niche_def` (resp. `gc`).
## 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`.
+In general, any type of temporal (resp. stratigraphic) data can be transformed. For this, you need define a `S3` class for the object you would like to transform (if it does not already have a S3 class), and then implement the transformation of this class 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 (CONTRIBUTING.md file) of the `admtools` package for details on how to contribute code.
+
+## References {#References}
-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.
+- Burgess, Peter. 2013. "CarboCAT: A cellular automata model of heterogeneous carbonate strata." Computers & Geosciences. .
-## References
+- Burgess, Peter. 2023. "CarboCATLite v1.0.1." Zenodo.
-\* 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.
+- 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.
diff --git a/vignettes/event_data.Rmd b/vignettes/event_data.Rmd
index 4a7a3b0..502728a 100644
--- a/vignettes/event_data.Rmd
+++ b/vignettes/event_data.Rmd
@@ -14,10 +14,9 @@ knitr::opts_chunk$set(
)
```
-
## Introduction
-This vignette provides an introduction to stratigraphic paleobiology applied to "_event type data_" such as location and age of specimens, and first/last occurrences of taxa.
+This vignette provides an introduction to stratigraphic paleobiology applied to "*event type data*" such as location and age of specimens, and first/last occurrences of taxa.
First, let's load the required packages:
@@ -28,50 +27,56 @@ library(admtools)
## Modeling fossil occurrences and first/last fossil occurrences
-Location and age of fossil specimens and first/last occurrences of taxa can be considered distinct events that occur at specific points in time or stratigraphic positions. Up to a few notable exceptions, the same principles apply to this "_event type data_", so we treat them together.
+Location and age of fossil specimens and first/last occurrences of taxa can be considered distinct events that occur at specific points in time or stratigraphic positions. Up to a few notable exceptions, the same principles apply to this "*event type data*", so we treat them together.
-We model events as _point processes_. In the `StratPal` package, we provide two functions to simulate point processes:
+We model events as *point processes*. In the `StratPal` package, we provide two functions to simulate point processes:
-* `p3`, short for Poisson point process, to simulate events that occur at a constant rate
+- `p3`, short for Poisson point process, to simulate events that occur at a constant rate
-* `p3_var_rate` for variable rate Poisson point processes to simulate events that occur at variable rates (e.g. extinction events etc.)
+- `p3_var_rate` for variable rate Poisson point processes to simulate events that occur at variable rates (e.g. extinction events)
-The usage of `p3` is straightforward.
+The usage of `p3` is straightforward:
```{r}
-# simulate fossil occurrences over one Myr with an average of 8 occurrences per Myr
+# simulate fossil occurrences over one Myr with an average of 15 occurrences per Myr
p3(rate = 15, from = 0, to = 1) |>
- hist(main = "Fossil occurrences",
+ hist(main = "Fossil abundance",
xlab = "Time [Myr]",
- ylab = "# Fossils")
+ ylab = "# Specimens")
```
+Note that due to to random fluctuations, the histogram can differ a lot from the constant rate that we simulated. This effect is more pronounced for low rates.
+
For `p3_var_rate`, you can pass either a function that specifies the rate to `x`:
```{r}
# return 100 occurrences by setting n parameter
-p3_var_rate(x = sin, from = 0, to = 9, n = 100) |>
- hist(xlab = "Time [Myr]", main = "Fossil occurrences")
# note that negative rates (where sin < 0) are ignored
+p3_var_rate(x = sin, from = 0, to = 9, n = 100) |>
+ hist(xlab = "Time [Myr]",
+ ylab= "# Specimens",
+ main = "Fossil abundance")
```
-Alternatively, you can pass it two vectors `x` and `y`. The rate is then determined by linear interpolation between these vectors. This is equivalent to using the function `approxfun(x,y, rule = 2)` as input for `x`
+Alternatively, you can pass it two vectors `x` and `y`. The rate is then determined by linear interpolation between these vectors. This is equivalent to using the function `approxfun(x,y, rule = 2)` as input for `x`.
```{r}
-# decline in fossil occurrences from 50 to 0 over 1 Myr
+# decline in last occurrences from 50 to 0 over 1 Myr
p3_var_rate(x = c(0,1), y = c(50, 0), from = 0, to = 1, f_max = 50) |>
- hist(xlab = "Time [Myr]", main = "Fossil occurrences")
+ hist(xlab = "Time [Myr]",
+ main = "Last occurrences",
+ ylab = "# Last occurrences")
```
-Keep in mind that it is only a matter of interpretation whether the "events" produced by `p3` and `p3_var_rate` are age/timing or locations of specimens from one taxon or of first/last occurrences of multiple taxa. The conceptual framework holds for all event type data.
+Keep in mind that it is only a matter of interpretation whether the "events" produced by `p3` and `p3_var_rate` are age/timing or locations of specimens from one taxon or of first/last occurrences of multiple taxa. The conceptual framework holds for any type of event type.
-**Task:** Model an "extinction event" with a background rate of 3 species disappearing per Myr, and a peak in extinction rate where the rate is 10-fold increased over an interval of 200 kyr. How distinct is the peak in the histogram, and how much do the results differe due to randomness when you run the analysis multiple times?
+**Task:** Model an "extinction event" with a background rate of 3 species disappearing per Myr, and a peak in extinction rate where the rate is 10-fold increased over an interval of 200 kyr. How distinct is the peak in the histogram, and how much do the results differ due to randomness when you run the simulation multiple times?
## Age-depth models
-Throughout this examples, we the age-depth models from 2 km and 12 km offshore in scenario A.
+Throughout this examples, we use the age-depth models from 2 km and 12 km offshore in scenario A from [Hohmann et al. (2024)](#References). See `?scenarioA` and `vignette("StratPal")` for details on scenario A.
-```{r, figures-side, fig.show="hold", out.width="50%", fig.align='left'}
+```{r, figures-side, fig.show="hold", out.width="50%"}
adm_2km = tp_to_adm(t = scenarioA$t_myr, # 2 km from shore
h = scenarioA$h_m[,"2km"],
T_unit = "Myr",
@@ -96,20 +101,20 @@ L_axis_lab()
title("Age-depth model 12 km from shore")
```
-**Task:** Plot the sea level curve for scenarion A, and the water depth curves at the two locations. How is this connected to the age-depth models, specifically to the presence and duration of gaps and low/high sedimentation rates? You can use `plot_sed_rate_t`, `get_hiat_list` and `get_hiat_duration` from `admtools`.
+**Task:** Plot the sea level curve for scenario A, and the water depth curves at the two locations. How is this connected to the age-depth models, specifically to the presence and duration of gaps and low/high sedimentation rates? You can use `plot_sed_rate_t`, `get_hiat_list` and `get_hiat_duration` from `admtools`.
-## Stratigraphic distortions of "event" type data
+## Stratigraphic distortions of "event" type data {#Stratigraphic-distortions}
-Any event type of data can directly be transformed using an age-depth model via `time_to_strat`. Depending on how the events are interpreted, we need to choose different options.
+Any event type data can directly be transformed using an age-depth model via `time_to_strat`. Depending on how the events are interpreted, we need to choose different options.
### Fossil specimens
Assuming the events are fossil specimens of one taxon, we can examine where specimens appear in the stratigraphic column:
```{r}
-p3(rate = 200, from = min_time(adm_2km), to = max_time(adm_2km)) |>
- time_to_strat(adm_2km, destructive = TRUE) |>
- hist(xlab = "Stratigraphic height [m]",
+p3(rate = 200, from = min_time(adm_2km), to = max_time(adm_2km)) |> # constant rate in time domain
+ time_to_strat(adm_2km, destructive = TRUE) |> # transform into depth domain
+ hist(xlab = "Stratigraphic height [m]", # plot
main = "Fossil abundance 2 km offshore",
ylab = "# Fossils",
breaks = seq(from = min_height(adm_2km), to = max_height(adm_2km), length.out = 20))
@@ -120,29 +125,32 @@ Note the option `destructive = TRUE` in `time_to_strat`. This is because we assu
Assuming the same rate of fossil specimens, we get a very different pattern 12 km from shore:
```{r}
-p3(rate = 200, from = min_time(adm_12km), to = max_time(adm_12km)) |>
- time_to_strat(adm_12km, destructive = TRUE) |>
- hist(xlab = "Stratigraphic height [m]",
+p3(rate = 200, from = min_time(adm_12km), to = max_time(adm_12km)) |> # same rate as 2 km from shore
+ time_to_strat(adm_12km, destructive = TRUE) |> # use different adm for transformation
+ hist(xlab = "Stratigraphic height [m]", # plot histogram
main = "Fossil abundance 12 km offshore",
ylab = "# Fossils",
breaks = seq(from = min_height(adm_12km), to = max_height(adm_12km), length.out = 20))
```
-Not only is the section much shorter (not even 20 m), but there are 2 distinct peaks at the bottom of the section and around 13 m. Comparing with the age-depth model shows that these locations correspond to times of extreme stratigraphic condensation (low sedimentation rates, indicated by a near-horizontal age-depth model). As a result, more time is represented per stratigraphic increment, leading to accumulations of fossil occurrences (Hohmann 2021).
+Not only is the section much shorter (not even 20 m), but there are 2 distinct peaks at the bottom of the section and around 13 m. Comparing with the age-depth model shows that these locations correspond to times of extreme stratigraphic condensation (low sedimentation rates, indicated by a near-horizontal age-depth model). As a result, more time is represented per stratigraphic increment, leading to fossil accumulations ([Hohmann 2021](#References)).
+
+**Task:** How does fossil abundance change along an onshore-offshore gradient when you assume the same rate of fossil occurrences in the time domain? Are there intervals where you can not recover spikes in fossil abundance, i.e. an increase in rate? How are these intervals linked to the external drivers of the platform architecture such as sea level and water depth?
### First/Last occurrences of common taxa
-In the fossil record, we can not directly observe origination and extinction rates of taxa, but only their first and last occurrences. This type of event data can also directly be transformed using `time_to_strat` when taxa are very common. This assumption will be explored further below.
+In the fossil record, we can not directly observe origination and extinction rates of taxa, but only their first and last occurrences. This type of event data can directly be transformed using `time_to_strat` when we assume taxa are very common. This assumption will be explored further below.
```{r}
-p3(rate = 200, from = min_time(adm_2km), to = max_time(adm_2km)) |>
- time_to_strat(adm_2km, destructive = FALSE) |>
- hist(xlab = "Stratigraphic height [m]",
+p3(rate = 200, from = min_time(adm_2km), to = max_time(adm_2km)) |> # constant rate of last occ
+ time_to_strat(adm_2km, destructive = FALSE) |> # non-destructive transformation!
+ hist(xlab = "Stratigraphic height [m]", # plot histogram
main = "Last occurrences 2 km offshore",
+ ylab = "# Last occurrences",
breaks = seq(from = min_height(adm_2km), to = max_height(adm_2km), length.out = 20))
```
-This would be the pattern of last occurrences in the section 2 km offshore when the rate of last occurrences is constant! The peaks coincide with the hiatuses, as the last occurrences of the species that disappear during the hiatus cluster at the hiatus surface. Note here that we used the option `destructive = FALSE` for `time_to_strat` to make sure the last occurrences appear at the location of the hiatus. This is why the assumption of abundant taxa is important: If taxa are rare, their last occurrence does not match their actual time of extinction (Signor-Lipps effect), and the last occurrences can be found further below the hiatus.
+This would be the pattern of last occurrences in the section 2 km offshore when the rate of last occurrences is constant in the time domain! The peaks coincide with the hiatuses, as the last occurrences of the species that disappear during the hiatus cluster at the hiatus surface ([Holland and Patzkowsky 2015](#References)). Note here that we used the option `destructive = FALSE` for `time_to_strat` to make sure last occurrences that happen during a hiatus appear at the hiatus surface. This is why the assumption of abundant taxa is important: If taxa are rare, their last occurrence does not match their actual time of extinction (Signor-Lipps effect, [Signor and Lipps 1982](#References)), and their last occurrences would be found below the hiatus surface. We explore the role of abundance on the location of last occurrences further in the section on [range offset](#Range-offset).
The patter 12 km offshore is very different:
@@ -151,126 +159,150 @@ p3(rate = 200, from = min_time(adm_12km), to = max_time(adm_12km)) |>
time_to_strat(adm_12km, destructive = TRUE) |>
hist(xlab = "Stratigraphic height [m]",
main = "Last occurrences 12 km offshore",
+ ylab = "# Last occurrences",
breaks = seq(from = min_height(adm_12km), to = max_height(adm_12km), length.out = 20))
```
-It almost perfectly matches the pattern of the fossil occurrences. This is because the offshore section is much more dominated by condensation (which influences both fossil occurrences and first/last occurrences equally), while the onshore section is dominated by hiatuses (which affect first/last occurrences). This shows it is important to understand how "event" type data is affected by hiatuses: is it preserved, or is it destroyed?
+It almost perfectly matches the pattern of fossil abundance. This is because the offshore section is much more dominated by condensation (which influences both fossil abundance and first/last occurrences equally), while the onshore section is dominated by hiatuses (which affect first/last occurrences). This shows that it is important to understand how "event" type data is affected by hiatuses: is it preserved, or is it destroyed?
+
+Naturally, the same modeling pipeline is valid for first occurrences of common taxa.
-### Range offset
+**Task:** How does the pattern of last occurrences change across an onshore-offshore gradient? If you observed a spike in last occurrences in one of the section, for what stratigraphic information do you need to look to exclude the possibility that the spike is a stratigraphic artifact?
+
+### Range offset {#Range-offset}
+
+Here we explore the idea of *range offset*, which is the offset between the last occurrence of a taxon and its actual time of extinction ([Holland and Patzkowsky 2002](#References)).
+
+First, we model all occurrences of the taxon up to its extinction in the time domain, and save it in a variable:
```{r}
t_ext = 1.5 # time of "true" extinction
-r = 30 # rate of fossil occurrences
+r = 30 # rate of fossil occurrences
# simulate rate fossil occurrences of taxon
f_occ = p3(r, from = min_time(adm_2km), to = t_ext)
-hist(f_occ)
+hist(f_occ,
+ xlab = "Time [Myr]",
+ ylab = "# Fossils",
+ main = "Fossil abundance",
+ breaks = seq(from = min_time(adm_2km), to = max_time(adm_2km), length.out = 20))
+```
+
+Then we can determine the stratigraphic position of the highest preserved fossil:
+
+```{r}
+highest_occ = f_occ |> # take fossil occ. in time domain
+ time_to_strat(adm_2km, destructive = TRUE) |> # transform into stratigraphic domain, destroying fossils that coincide with hiatuses
+ max(na.rm = TRUE) # find highest preserved fossil (destroyed fossils are NA)
+```
-# stratigraphically highest fossil found
-highest_occ = f_occ |>
- time_to_strat(adm_2km, destructive = TRUE) |> # destroy fossil occurrences
- max(na.rm = TRUE)
+Now we can determine two types of range offset. The first is the distance in meters between the position of the last preserved fossil and the position where the taxon actually goes extinct:
-# stratigraphic position of "true" extinction
-h_true_ext = t_ext |>
+```{r}
+h_true_ext = t_ext |> # stratigraphic position of "true" extinction
time_to_strat(adm_2km, destructive = FALSE)
-# distance between last occurence of taxon and location where they actually go extinct
+# distance between last occurrence of taxon and location where they actually go extinct
strat_range_offset_m = h_true_ext - highest_occ
+```
-# time when last preserved fossil lived
-t_last_occ = highest_occ |>
- strat_to_time(adm_2km)
+In this example, the stratigraphic range offset is `r signif(strat_range_offset_m, digits = 3)` meters.
+
+The second type of range offset is the difference in time between the age of the last occurrence and the actual time of extinction:
+```{r}
+t_last_occ = highest_occ |> # time when last preserved fossil lived
+ strat_to_time(adm_2km)
# time offset between true extinction and time when last fossil lived.
time_range_offset_myr = t_ext - t_last_occ
```
+In this example, the temporal range offset is `r signif(time_range_offset_myr, digits = 3)` Myr years.
+**Task:** How does the (temporal and stratigraphic) range offset change with the rate parameter? For which rate values are taxa abundant enough to make the Signor-Lipps effect negligible (i.e. the stratigraphic range offset is sufficiently small)? How do the results vary with (1) different times of "true" extinction within a fixed section and (2) across the onshore-offshore gradient?
-```{r}
-dist = scenarioA$dist_from_shore[2]
+## Niche modeling
-# define age-depth model and attach units to it
-adm = tp_to_adm(t = scenarioA$t_myr,
- h = scenarioA$h_m[,dist],
- T_unit = "Myr",
- L_unit = "m")
-# plot age-depth model
-plot(adm,
- lwd_acc = 2, # plot thicker lines for accumulative intervals (lwd = line width)
- lty_destr = 0) # don't plot destructive intervals (lty = line type)
-T_axis_lab() # add time axis label
-L_axis_lab() # add length axis label
-```
+So far, we have only examined the effects of stratigraphic distortions (i.e., gaps and changes in sedimentation rates) on event type data. Here, we increase the complexity of the pipeline by incorporating niche modeling. For this, we need two components:
-```{r}
-p3(rate = 100, from = min_time(adm), to = max_time(adm)) |> # simulate fossil occurrences in time
- get_height(adm, t = _, destructive = FALSE) |> # transfrom into depth domain
- hist(breaks = seq(min_height(adm), to = max_height(adm), length.out = 20),
- main = "Fossil occurrences",
- xlab = "Height [m]")
-```
+- a niche definition that describes how "comfortable" a taxon feels along an (environmental or eclogical) gradient
+- information on how that gradient changes with time
-## Influence of niches
+In this example we focus on *water depth* as gradient. The general modeling framework works for any type of gradient, see the "advanced functionality" vignette for details.
-We increase the complexity of the pipeline by incorporating niche modeling. For this we need two components:
+First, we define a function that determines how water depth (our gradient) changes with time:
-* a niche definition that describes how "comfortable" a species feels along an environmental gradient
+```{r}
+t = scenarioA$t_myr # time steps of the model
+wd = scenarioA$wd_m[,"2km"] # water depth 2 km offshore at model time steps
+gc = approxfun(t, wd) # define function that defines how the gradient changes with time (gc = *G*radient *C*hange)
+
+plot(t, gc(t),
+ type = "l",
+ xlab = "Time [Myr]",
+ ylab = "Water depth [m]",
+ main = "Water depth 2 km offshore")
+```
-* information on how that gradient changes with time
+Next, we define the niche. This is done using a function that takes water depth (our gradient) as input, and returns numbers between 0 and 1. Here a 1 reflects that the taxon is fully within its preferred niche, and a 0 reflects that it is completely outside of its preferred niche. Values between 0 and 1 are interpreted as probabilities that the taxon is observed locally.
-In this example we focus on _water depth_ as gradient as it is included in the model outputs. The general modeling framework works for any type of gradient.
+For our example, we use the helper function `cnd` (capped normal distribution) to define a niche. It is a probability density function of the normal distribution, multiplied by the factor `inc` and capped at 1. This is identical to the *probability of collection* model by [Holland and Patzkowsky (2002)](#References). See the "advanced functionality" vignette for details how to construct other niche definitions.
-First, we define the niche. We do this using the helper function `cnd` (capped normal distribution). It is a probability density function of the normal distribution, multiplied by the factor `inc` and capped at 1. A 1 reflects that the taxon is fully within its prefferred environment, and a 0 reflects that is is outside of its environment.
```{r}
-my_niche = cnd(mean = 10, # preffered water depth
- sd = 5, # tolerance to water depth fluctiations
- inc = 15, # multiplier
- cut_neg = TRUE) # cut off negative avalues - the species does not surviva at land
+my_niche = cnd(mean = 10, # preferred water depth
+ sd = 5, # depth tolerance
+ inc = 15, # multiplier
+ cut_neg = TRUE) # cut off negative values - the taxon does not survive at land
x = seq(-2, 30, by = 0.1)
-plot(x, my_niche(x), type = "l", xlab = "water depth", ylab = "preference")
+plot(x, my_niche(x),
+ type = "l",
+ xlab = "Water depth [m]",
+ ylab = "Preference",
+ main = "Water depth preference")
```
-To get the gradient change with time, we extract water depth 2 km offshore from the example data:
+The taxon we are modeling prefers to live in 10 m water depth (`mean = 10`), is moderately tolerant to fluctuations in water depth (`sd = 5`) and does not survive on land (`cut_neg = TRUE`).
-```{r}
-t = scenarioA$t_myr
-wd = scenarioA$wd_m[,"2km"]
-gc = approxfun(t, wd) # define gradient change
-plot(t, gc(t), type = "l", xlab = "Time [Myr]", ylab = "Water depth [m]")
-
-```
-
-Now we can combine these two ingredients to examine how the preffered niche changes the fossil abundance in the time domain.
+With the niche and the change in gradient defined, we can use `apply_niche_pref` for niche modeling.
```{r}
-p3(rate = 300, from = min_time(adm_2km), to = max_time(adm_2km)) |>
- apply_niche_pref(niche_def = my_niche, gc = gc) |>
- hist(xlab = "Time [Myr]", main = "Fossil abundance")
+p3(rate = 300, from = min_time(adm_2km), to = max_time(adm_2km)) |> # model occurrences of taxon based on a constant rate
+ apply_niche_pref(niche_def = my_niche, gc = gc) |> # apply the niche preferenc of the taxon as defined above
+ hist(xlab = "Time [Myr]",
+ main = "Fossil abundance",
+ ylab = "# Fossils",
+ breaks = seq(from = min_time(adm_2km), to = max_time(adm_2km), length.out = 40))
```
Here you can see that the taxon only occurs at times where the water depth is favorable for it. Around 2 Myr, the water depth is too shallow, and the fossil abundance is reduced.
-With the basic niche modeling done, we can now combine it with the stratigrahic distortion by age-depth models to get an idea how niche preferences change the fossil abundance in the section.
+With the basic niche modeling done, we can incorporate it into our modeling pipeline to examine the joint effects of stratigraphic distortion by age-depth models and niche preferences.
```{r}
-p3(rate = 300, from = min_time(adm_2km), to = max_time(adm_2km)) |>
- apply_niche_pref(niche_def = my_niche, gc = gc) |>
- time_to_strat(adm_2km, destructive = TRUE) |>
- hist(xlab = "Stratigraphic height [m]",
- main = "Fossil abundance")
+p3(rate = 300, from = min_time(adm_2km), to = max_time(adm_2km)) |> # model occurrences based on constant rate
+ apply_niche_pref(niche_def = my_niche, gc = gc) |> # apply niche preference
+ time_to_strat(adm_2km, destructive = TRUE) |> # transform into strat. domain, destroy fossils that coincide with hiatuses
+ hist(xlab = "Stratigraphic height [m]", # plot results
+ main = "Fossil abundance 2 km from shore",
+ ylab = "# Fossils",
+ breaks = seq(from = min_height(adm_2km), to = max_height(adm_2km), length.out = 40))
```
-We can see a general decrease in fossil abundance, but no longer gaps as were visible in the time domain. The latter is because the gaps in fossil abundance in the time domain coincide with hiatuses, and will thus not be visible in the stratigraphic domain. The general decrease in fossil abundance is because as the carbonate platform ages, the water depth on the platform is decreasing:
+In most cases, you can see a general decrease in fossil abundance, but in contrast to the time domain, no prolonged intervals without fossil are visible. This latter is because the gaps in fossil abundance in the time domain coincide with hiatuses, and will thus not be visible in the stratigraphic domain. The general decrease in fossil abundance is because as the carbonate platform ages, the water depth on top of the platform decreases:
```{r}
-list("t" = t, "y" = wd) |>
- time_to_strat(adm_2km) |>
- plot(orientation = "lr", type = "l", xlab = "Stratigraphic position [m]", ylab = "Water depth [m]")
+list("t" = t, "y" = wd) |> # create list with time - water depth information
+ time_to_strat(adm_2km) |> # transform into the strat. domain
+ plot(orientation = "lr", # plot water depth information in the stratigraphic domain
+ type = "l",
+ xlab = "Stratigraphic position [m]",
+ ylab = "Water depth [m]",
+ main = "Water depth in section 2 km from shore")
```
+**Task:** How does the effect of niche preference change across the onshore-offshore gradient? There is a dependency between stratigraphic effects (gaps, changes in sedimentation rate) and water depth/eustatic sea level. Does linking ecology to water depth strengthen or weaken the effects of stratigraphic distortions examined in section ["stratigraphic distortions of event type data"](#Stratigraphic-distortions)?
+
## Further reading
Go to
@@ -279,7 +311,7 @@ Go to
vignette("phenotypic_evolution")
```
-for details on how to model stratigraphic paleobiology of phenotypic evolution, or explore the vignette online under ([mindthegap-erc.github.io/StratPal/articles/phenotypic_evolution](https://mindthegap-erc.github.io/StratPal/articles/phenotypic_evolution.html)).
+for details on how to model stratigraphic paleobiology of phenotypic evolution, or explore the vignette online under [mindthegap-erc.github.io/StratPal/articles/phenotypic_evolution](https://mindthegap-erc.github.io/StratPal/articles/phenotypic_evolution.html).
See also
@@ -287,8 +319,16 @@ See also
vignette("advanced_functionality")
```
-for details on how to expand on the modeling pipelines described here, or explore the vignette online under ([mindthegap-erc.github.io/StratPal/articles/advanced_functionality](https://mindthegap-erc.github.io/StratPal/articles/advanced_functionality.html)).
+for details on how to expand on the modeling pipelines described here, or explore the vignette online under [mindthegap-erc.github.io/StratPal/articles/advanced_functionality](https://mindthegap-erc.github.io/StratPal/articles/advanced_functionality.html).
+
+## References {#References}
+
+* Hohmann, Niklas. 2021. "Incorporating information on varying sedimentation rates into paleontological analyses." PALAIOS. https://doi.org/10.2110/palo.2020.038.
+
+* 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.
+
+* 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
-## Reference
+* Holland, Steven M. and Patzkowsky, Mark E. 2015. "The stratigraphy of mass extinction." Palaeontology. https://doi.org/10.1111/pala.12188.
-* NIKLAS HOHMANN; INCORPORATING INFORMATION ON VARYING SEDIMENTATION RATES INTO PALEONTOLOGICAL ANALYSES. PALAIOS 2021;; 36 (2): 53–67. doi: https://doi.org/10.2110/palo.2020.038
+* Signor, Philip W.; Lipps, Jere H. 1982. "Sampling bias, gradual extinction patterns and catastrophes in the fossil record." Geological Society of America Special Paper 190.
diff --git a/vignettes/phenotypic_evolution.Rmd b/vignettes/phenotypic_evolution.Rmd
index 095c2ef..6fccdd7 100644
--- a/vignettes/phenotypic_evolution.Rmd
+++ b/vignettes/phenotypic_evolution.Rmd
@@ -57,7 +57,7 @@ We are interested in how phenotypic evolution is preserved in the stratigraphic
### Age-depth models
-As an example, we compare the preservation of trait evolution 2 km and 12 km offshore in the carbonate platform in scenario A from [Hohmann et al. 2024](#References). See `?scenarioA` and `vignette("StratPal")` for details on scenario A.
+As an example, we compare the preservation of trait evolution 2 km and 12 km offshore in the carbonate platform in scenario A from [Hohmann et al. (2024)](#References). See `?scenarioA` and `vignette("StratPal")` for details on scenario A.
First define the age-depth models:
From 3b6581e0949a2030fc1291a410e490a9b35cc00f Mon Sep 17 00:00:00 2001
From: Niklas Hohmann <67792281+NiklasHohmann@users.noreply.github.com>
Date: Sat, 20 Jul 2024 17:24:28 +0200
Subject: [PATCH 5/5] expand documentation
---
R/StratPal-package.R | 6 ++++++
R/apply_niche_pref.R | 14 +++++++++-----
R/apply_taphonomy.R | 13 +++++++------
R/cnd.R | 10 ++++++++--
R/cnd_niche.R | 17 -----------------
R/ornstein_uhlenbeck.R | 12 ++++++------
R/p3.R | 10 +++++++++-
R/p3_var_rate.R | 9 +++++++--
R/random_walk.R | 12 ++++++------
R/rej_samp.R | 4 ++--
R/scenarioA.R | 19 +++++++++++++++----
R/stasis.R | 9 +++++----
R/thin.R | 9 +++++----
man/StratPal-package.Rd | 24 ++++++++++++++++++++++++
man/apply_niche_pref.Rd | 14 +++++++++-----
man/apply_taphonomy.Rd | 12 ++++++------
man/cnd.Rd | 14 ++++++++++++--
man/ornstein_uhlenbeck.Rd | 12 ++++++------
man/p3.Rd | 12 +++++++++++-
man/p3_var_rate.Rd | 11 +++++++++--
man/random_walk.Rd | 12 ++++++------
man/rej_samp.Rd | 6 ++++--
man/scenarioA.Rd | 16 +++++++++++++---
man/stasis.Rd | 11 +++++++----
man/thin.Rd | 9 +++++----
25 files changed, 197 insertions(+), 100 deletions(-)
create mode 100644 R/StratPal-package.R
delete mode 100644 R/cnd_niche.R
create mode 100644 man/StratPal-package.Rd
diff --git a/R/StratPal-package.R b/R/StratPal-package.R
new file mode 100644
index 0000000..a65cf64
--- /dev/null
+++ b/R/StratPal-package.R
@@ -0,0 +1,6 @@
+#' @keywords internal
+"_PACKAGE"
+
+## usethis namespace: start
+## usethis namespace: end
+NULL
diff --git a/R/apply_niche_pref.R b/R/apply_niche_pref.R
index d0a2eab..7d77233 100644
--- a/R/apply_niche_pref.R
+++ b/R/apply_niche_pref.R
@@ -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{
@@ -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))
diff --git a/R/apply_taphonomy.R b/R/apply_taphonomy.R
index 3d500df..a606e8a 100644
--- a/R/apply_taphonomy.R
+++ b/R/apply_taphonomy.R
@@ -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)
diff --git a/R/cnd.R b/R/cnd.R
index a9240e7..c88223a 100644
--- a/R/cnd.R
+++ b/R/cnd.R
@@ -13,7 +13,7 @@ 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
@@ -21,11 +21,17 @@ cnd = function(mean, sd, inc = 1, cut_neg = TRUE){
#' # 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)))
diff --git a/R/cnd_niche.R b/R/cnd_niche.R
deleted file mode 100644
index 73d1bdb..0000000
--- a/R/cnd_niche.R
+++ /dev/null
@@ -1,17 +0,0 @@
-#' cnd = function(mean, sd, inc = 1){
-#' #' @export
-#' #'
-#' #' @title capped normal distribution pdf
-#' #'
-#' #' @param mean mean of normal distribution
-#' #' @param sd standard deviation
-#' #' @param inc scalar, factor by which the pdf in multiplied
-#' #'
-#' #' @returns a function
-#' #'
-#' #' @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)
-#'
-#' f = function(x) pmin(inc * dnorm(x, mean, sd), rep(1, length(x)))
-#' return(f)
-#' }
diff --git a/R/ornstein_uhlenbeck.R b/R/ornstein_uhlenbeck.R
index 1d40c9b..51ceada 100644
--- a/R/ornstein_uhlenbeck.R
+++ b/R/ornstein_uhlenbeck.R
@@ -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{
diff --git a/R/p3.R b/R/p3.R
index 46b4e32..607c0e1 100644
--- a/R/p3.R
+++ b/R/p3.R
@@ -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\".")}
diff --git a/R/p3_var_rate.R b/R/p3_var_rate.R
index 73a393e..dba911e 100644
--- a/R/p3_var_rate.R
+++ b/R/p3_var_rate.R
@@ -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\".")
diff --git a/R/random_walk.R b/R/random_walk.R
index fa47730..decf80c 100644
--- a/R/random_walk.R
+++ b/R/random_walk.R
@@ -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{
diff --git a/R/rej_samp.R b/R/rej_samp.R
index 7913f62..b348ce2 100644
--- a/R/rej_samp.R
+++ b/R/rej_samp.R
@@ -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
@@ -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.")}
diff --git a/R/scenarioA.R b/R/scenarioA.R
index a64c6b8..70008bb 100644
--- a/R/scenarioA.R
+++ b/R/scenarioA.R
@@ -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:
#'
@@ -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. .
+ #' * Burgess, Peter. 2023. "CarboCATLite v1.0.1." Zenodo.
+ #' * 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"
diff --git a/R/stasis.R b/R/stasis.R
index 36ca1cc..dfa8783 100644
--- a/R/stasis.R
+++ b/R/stasis.R
@@ -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
diff --git a/R/thin.R b/R/thin.R
index 165c494..995b280 100644
--- a/R/thin.R
+++ b/R/thin.R
@@ -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
diff --git a/man/StratPal-package.Rd b/man/StratPal-package.Rd
new file mode 100644
index 0000000..5eedace
--- /dev/null
+++ b/man/StratPal-package.Rd
@@ -0,0 +1,24 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/StratPal-package.R
+\docType{package}
+\name{StratPal-package}
+\alias{StratPal}
+\alias{StratPal-package}
+\title{StratPal: Stratigraphic Paleobiology Modeling Pipelines}
+\description{
+Models for stratigraphic paleobiology, including niche modeling and age-depth modeling.
+}
+\seealso{
+Useful links:
+\itemize{
+ \item \url{https://mindthegap-erc.github.io/StratPal/ }
+ \item \url{https://github.com/MindTheGap-ERC/StratPal}
+ \item Report bugs at \url{https://github.com/MindTheGap-ERC/StratPal/issues}
+}
+
+}
+\author{
+\strong{Maintainer}: Niklas Hohmann \email{N.H.Hohmann@uu.nl} (\href{https://orcid.org/0000-0003-1559-1838}{ORCID})
+
+}
+\keyword{internal}
diff --git a/man/apply_niche_pref.Rd b/man/apply_niche_pref.Rd
index 3f74b66..2a24a43 100644
--- a/man/apply_niche_pref.Rd
+++ b/man/apply_niche_pref.Rd
@@ -7,15 +7,15 @@
apply_niche_pref(x, niche_def, gc)
}
\arguments{
-\item{x}{events, e.g. fossil occurrences}
+\item{x}{events, e.g. times/ages of fossil occurrences or their stratigraphic position.}
-\item{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.}
+\item{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.}
\item{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 \code{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 \code{niche_def} and \code{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 \code{thin}.
+Combines the functions \code{niche_def} and \code{gc} ("gradient change") to determine how the taxons' preferred environment changes with time. This is done by composing \code{niche_def} and \code{gc}. The result is then used as a thinning on the events \code{x}.
}
\examples{
\dontrun{
@@ -42,9 +42,13 @@ Combines the functions niche_def and gc ("gradient change") to determine how the
# fossil occurrences with niche preference
hist(foss_occ_niche, xlab = "time")
+ # see also
+ vignette("event_data")
+ # for a longer example
+
}
}
\seealso{
-\code{\link[=apply_taphonomy]{apply_taphonomy()}} to model taphonomic effects based on the sampe principle, \code{\link[=thin]{thin()}} for the underlying mathematical procedure
+\code{\link[=apply_taphonomy]{apply_taphonomy()}} to model taphonomic effects based on the same principle, \code{\link[=thin]{thin()}} for the underlying mathematical procedure.
}
diff --git a/man/apply_taphonomy.Rd b/man/apply_taphonomy.Rd
index c2a63ec..369054a 100644
--- a/man/apply_taphonomy.Rd
+++ b/man/apply_taphonomy.Rd
@@ -2,20 +2,20 @@
% Please edit documentation in R/apply_taphonomy.R
\name{apply_taphonomy}
\alias{apply_taphonomy}
-\title{model taphonomy}
+\title{Model taphonomic effects}
\usage{
apply_taphonomy(x, pres_potential, ctc)
}
\arguments{
-\item{x}{numeric vector. event type data, e.g. fossil occurrences}
+\item{x}{events, e.g. times/ages of fossil occurrences or their stratigraphic position.}
-\item{pres_potential}{function. takes taphonomic conditions as input and returns the preservation potential (a number between 0 and 1)}
+\item{pres_potential}{function. Takes taphonomic conditions as input and returns the preservation potential (a number between 0 and 1)}
-\item{ctc}{function, change in taphonomic conditions (ctc) with time.}
+\item{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 \code{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 \code{thin}.
}
\seealso{
-\code{\link[=apply_niche_pref]{apply_niche_pref()}} for modeling niche preferences based on the same principle, \code{\link[=thin]{thin()}} for the underlying mathematical procedure
+\code{\link[=apply_niche_pref]{apply_niche_pref()}} for modeling niche preferences based on the same principle, \code{\link[=thin]{thin()}} for the underlying mathematical procedure.
}
diff --git a/man/cnd.Rd b/man/cnd.Rd
index 8af120d..3b9602b 100644
--- a/man/cnd.Rd
+++ b/man/cnd.Rd
@@ -19,11 +19,12 @@ cnd(mean, sd, inc = 1, cut_neg = TRUE)
a function
}
\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).
}
\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
@@ -31,7 +32,16 @@ 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
}
}
+\references{
+\itemize{
+\item 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
+}
+}
diff --git a/man/ornstein_uhlenbeck.Rd b/man/ornstein_uhlenbeck.Rd
index 0fb21d6..7e8dcb6 100644
--- a/man/ornstein_uhlenbeck.Rd
+++ b/man/ornstein_uhlenbeck.Rd
@@ -7,21 +7,21 @@
ornstein_uhlenbeck(t, mu = 0, theta = 1, sigma = 1, y0 = 0)
}
\arguments{
-\item{t}{times at which the process is simulated}
+\item{t}{times at which the process is simulated. Can be heterodistant}
-\item{mu}{numeric, long term mean}
+\item{mu}{scalar, long term mean}
-\item{theta}{numeric, mean reversion speed}
+\item{theta}{scalar, mean reversion speed}
\item{sigma}{positive scalar, strength of randomness}
-\item{y0}{scalar. initial value (value of process at the first entry of t)}
+\item{y0}{scalar, initial value (value of process at the first entry of t)}
}
\value{
-a list with two elements: \code{t} and \code{y}. \code{t} is a duplicate of the input \code{t}, \code{y} are the values of the OU process at these times. Outputs are of class \code{timelist} and can thus be plotted directly using \code{plot}, see \code{?plot.timelist}
+A list with two elements: \code{t} and \code{y}. \code{t} is a duplicate of the input \code{t}, \code{y} are the values of the OU process at these times. Output list is of class \code{timelist} and can thus be plotted directly using \code{plot}, see \code{?admtools::plot.timelist}
}
\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 \code{t}.
+Simulates an Ornstein-Uhlenbeck process using the Euler-Maruyama method. The process is simulated on a scale of \code{0.25 * min(diff(t))} and then interpolated to the values of \code{t}.
}
\examples{
\dontrun{
diff --git a/man/p3.Rd b/man/p3.Rd
index 0195db2..87bf87c 100644
--- a/man/p3.Rd
+++ b/man/p3.Rd
@@ -15,6 +15,9 @@ p3(rate, from, to, n = NULL)
\item{n}{integer of NULL (default). Number of events to return. If NULL, the number is random and determined by the rate parameter}
}
+\value{
+a numeric vector with timing/location of events.
+}
\description{
Simulates events in the interval \code{from} to \code{to} based on a Poisson point process with rate \code{rate}. If the parameter \code{n} is used, the number of fossils is conditioned to be \code{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
@@ -23,13 +26,20 @@ In the context of paleontology, these events can be interpreted as fossil occurr
\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{
+\code{\link[=p3_var_rate]{p3_var_rate()}} for the variable rate implementation
+}
diff --git a/man/p3_var_rate.Rd b/man/p3_var_rate.Rd
index 0bd25de..16555e1 100644
--- a/man/p3_var_rate.Rd
+++ b/man/p3_var_rate.Rd
@@ -27,16 +27,23 @@ In the context of paleontology, these events can be interpreted as fossil occurr
\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{
+\code{\link[=p3]{p3()}} for the constant rate implementation, \code{\link[=rej_samp]{rej_samp()}} for the underlying random number generation.
+}
diff --git a/man/random_walk.Rd b/man/random_walk.Rd
index 4bca43a..dc37d21 100644
--- a/man/random_walk.Rd
+++ b/man/random_walk.Rd
@@ -7,19 +7,19 @@
random_walk(t, sigma = 1, mu = 0, y0 = 0)
}
\arguments{
-\item{t}{numeric vector with strictly increasing elements. Times at which the random walk is evaluated}
+\item{t}{numeric vector with strictly increasing elements, can be heterodistant. Times at which the random walk is evaluated}
-\item{sigma}{variance parameter}
+\item{sigma}{positive scalar, variance parameter}
-\item{mu}{directionality parameter}
+\item{mu}{scalar, directionality parameter}
-\item{y0}{starting value, i.e. value of the random walk at the first entry of \code{t}}
+\item{y0}{scalar, starting value (value of the random walk at the first entry of \code{t})}
}
\value{
-a list with elements \code{t} and \code{y}. \code{t} is a duplicate of the input parameter and is the times at which the random walk is evaluated. \code{y} are the values of the random walk at said times
+A list with elements \code{t} and \code{y}. \code{t} is a duplicate of the input parameter and is the times at which the random walk is evaluated. \code{y} are the values of the random walk at said times. Output list is of class \code{timelist} and can thus be plotted directly using \code{plot}, see \code{?admtools::plot.timelist}
}
\description{
-simulates a random walk as a Brownian drift
+Simulates a (continuous time) random walk as a Brownian drift
}
\examples{
\dontrun{
diff --git a/man/rej_samp.Rd b/man/rej_samp.Rd
index 3495f22..8ba6312 100644
--- a/man/rej_samp.Rd
+++ b/man/rej_samp.Rd
@@ -18,7 +18,7 @@ rej_samp(f, x_min, x_max, n = 1L, f_max = 1)
\item{f_max}{maximum value of \code{f} in the interval from \code{x_min} to \code{x_max}. If f attains values larger than \code{f_max} a warning is throw, \code{f_max} is adjusted, and sampling is started again}
}
\description{
-rejection sampling from the (pseudo) pdf \code{f} in the interval between \code{x_min} and \code{x_max}. Returns \code{n} samples. Note that values of \code{f} below 0 are capped to zero
+Rejection sampling from the (pseudo) pdf \code{f} in the interval between \code{x_min} and \code{x_max}. Returns \code{n} samples. Note that values of \code{f} below 0 are capped to zero
}
\examples{
\dontrun{
@@ -26,5 +26,7 @@ f = sin
x = rej_samp(f, 0, 3*pi, n = 100)
hist(x) # note that no samples are drawn where sin is negative
}
-
+}
+\seealso{
+\code{\link[=p3_var_rate]{p3_var_rate()}} for the derived variable rate Poisson point process implementation.
}
diff --git a/man/scenarioA.Rd b/man/scenarioA.Rd
index 8a9300c..038777f 100644
--- a/man/scenarioA.Rd
+++ b/man/scenarioA.Rd
@@ -3,7 +3,7 @@
\docType{data}
\name{scenarioA}
\alias{scenarioA}
-\title{scenario A from Hohmann et al 2024}
+\title{example data, scenario A from Hohmann et al. (2024)}
\format{
A list with 6 elements:
\itemize{
@@ -12,15 +12,25 @@ A list with 6 elements:
\item \code{dist_from_shore} : character vector. Distance from shore in km of locations at which the observations were made
\item \code{h_m} : matrix of size length(t_myr) x length(dist_from_shore). Accumulated sediment height in m at examined locations
\item \code{wd_m}: matrix of size length(t_myr) x length(dist_from_shore). Water depth in m at examined locations
-\item \code{strat_col}: list with length(dist_from shore) elements. Represents a stratigraphic column. Each element is a list with two elements
+\item \code{strat_col}: list with length(dist_from shore) elements. Represents a stratigraphic column. Each element is a list with two elements:
+\itemize{
\item \code{bed_thickness_m}: numeric vector. Bed thickness in m
\item \code{facies_code} : integer vector. facies code of the bed
}
}
+}
\usage{
scenarioA
}
\description{
-scenario A from Hohmann et al 2024
+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)
+}
+\references{
+\itemize{
+\item Burgess, Peter. 2013. "CarboCAT: A cellular automata model of heterogeneous carbonate strata." Computers & Geosciences. \url{https://doi.org/10.1016/j.cageo.2011.08.026}.
+\item Burgess, Peter. 2023. "CarboCATLite v1.0.1." Zenodo. \url{https://doi.org/10.5281/zenodo.8402578}
+\item 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.
+\item 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 \href{https://creativecommons.org/licenses/by/4.0/}{CC-BY 4.0} license.
+}
}
\keyword{datasets}
diff --git a/man/stasis.Rd b/man/stasis.Rd
index ff0d74d..f169ad4 100644
--- a/man/stasis.Rd
+++ b/man/stasis.Rd
@@ -7,14 +7,17 @@
stasis(t, mean = 0, sd = 1)
}
\arguments{
-\item{t}{times at which the phenotype is determined}
+\item{t}{times at which the traits are determined}
-\item{mean}{mean value}
+\item{mean}{scalar, mean trait value}
-\item{sd}{standard deviation}
+\item{sd}{strictly positive scalar, standard deviation of traits}
+}
+\value{
+A list with two elements: \code{t} and \code{y}. \code{t} is a duplicate of the input \code{t}, \code{y} are the corresponding trait values. Output list is of class \code{timelist} and can thus be plotted directly using \code{plot}, see \code{?admtools::plot.timelist}
}
\description{
-simulates stasis as independent, normally distributed random variables with mean \code{mean} and standard deviation \code{sd}
+Simulates stasis as independent, normally distributed random variables with mean \code{mean} and standard deviation \code{sd}
}
\examples{
\dontrun{
diff --git a/man/thin.Rd b/man/thin.Rd
index 060232b..f3673e5 100644
--- a/man/thin.Rd
+++ b/man/thin.Rd
@@ -12,17 +12,18 @@ thin(x, thin)
\item{thin}{a function used for thinning}
}
\description{
-Thins a vector of events using the function thin, meaning the probability that the ith event in x is preserved is given by \emph{thin (x(i))}. Values of
-\code{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, \code{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 \emph{thin(x(i))}. Values of
+\code{thin} below 0 and above 1 are ignored.
+Is used to model niche preferences in \code{apply_niche_pref} and taphonomic effects in \code{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{