---
title: "Cleaning fossil records"
description: >
  The age tests and PyRate export for Paleobiology Database records.
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Cleaning fossil records}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>", fig.width = 6.5, fig.height = 4,
                      dev = "svglite", fig.ext = "svg", dev.args = list(pointsize = 14))
has_cc <- requireNamespace("CoordinateCleaner", quietly = TRUE)
```

Fossil occurrences, such as those downloaded from the Paleobiology Database (PBDB), can carry the same coordinate errors as modern records, and they have one more field to check: the age of the record, given as a range in million years (Ma). This article flags a small fossil table with the coordinate tests and the age test `age_equal`, filters records by the width of their age range, tests a larger table for outliers within a taxon, and exports the kept records to PyRate. [Get started](OccurrenceCleaner.html) covers installation; the comparison with [CoordinateCleaner](https://ropensci.github.io/CoordinateCleaner/) at the end also needs that package installed.

1. **Flag** a fossil table with `flag_occurrences()`.

2. **Filter** out flagged records, using `keep_unflagged()`.

3. **Filter** by the width of the age range, counting each group with `tally_flags()`.

4. **Test** each taxon for outliers in age and space.

5. **Export** the kept records with `export_pyrate()`.

## A small fossil table

PBDB uses `accepted_name` for taxon names, `lng` and `lat` for coordinates, and `cc` for two-letter country codes. The columns `min_ma` and `max_ma` give the youngest and oldest possible ages in millions of years before the present.

The example below contains nine synthetic records, each built for one case; the cases are explained beside the flags they produce. Records 1 and 9 are the controls: a minimum age below the maximum age, a span under 4 million years, and coordinates in Wyoming that match their country code `US`. Records 7 and 8 have the country code `DE`; their coordinates are missing in the table and set in the next chunk:

```{r data}
library(OccurrenceCleaner)

fos <- data.frame(
  occurrence_no = 1:9,
  accepted_name = c("Coryphodon", "Coryphodon", "Hyracotherium", "Hyracotherium",
                    "Hyracotherium", "Mesonyx", "Mesonyx", "Mesonyx", "Uintatherium"),
  lng = c(-108.52, -108.40, -108.71, -107.95, -108.10, 0, NA, NA, -109.03),
  lat = c(44.31, 44.05, 43.92, 44.60, 44.40, 0, NA, NA, 41.62),
  cc = c("US", "US", "US", "US", "US", "US", "DE", "DE", "US"),
  min_ma = c(53.0, 55.8, 56.0, 23.0, NA, 48.6, 47.8, 47.8, 46.2),
  max_ma = c(56.0, 55.8, 52.5, 66.0, 54.0, 50.2, 49.0, 49.0, 48.6))
```

Records 7 and 8 are placed relative to the centroid of Germany in the package's reference points, so their distance to it is known: record 7 sits on the point, record 8 lies 0.027 degrees of latitude (about 3 km) north of it. Their ages are ordered, with a span of 1.2 million years.

```{r centroid}
cen <- reference_points("centroid")
de <- cen[cen$key == "ne:DEU" & cen$method == "polygon centroid", ]
de[, c("name", "method", "lon", "lat")]
fos[7:8, "lng"] <- de$lon
fos[7:8, "lat"] <- de$lat + c(0, 0.027)
```

## Flagging

`occ_columns()` maps the PBDB names onto the fields the tests read. The coordinate tests are the record tests of [`clean_fossils()`](https://ropensci.github.io/CoordinateCleaner/reference/clean_fossils.html) in CoordinateCleaner; `age_equal` is its `agesequal` test, and the country check runs too because the table has country codes:

```{r flag}
cols <- occ_columns(lon = "lng", lat = "lat", species = "accepted_name",
                    country = "cc", id = "occurrence_no")
tests <- c("invalid", "zero", "equal", "centroid", "institution",
           "gbif_headquarters", "country_mismatch", "age_equal")
f <- flag_occurrences(fos, tests = tests, columns = cols)
f[, c("occurrence_no", "accepted_name", "min_ma", "max_ma", "flags")]
```

`age_equal` flags a record whose minimum age is equal to or greater than its maximum age. Record 2 has equal minimum and maximum ages and record 3 has them reversed, so both are flagged. A zero-width range can come from one age copied into both fields, and an inverted one from swapped fields, but the flag only records the comparison; such records need a look at their source before they are corrected or dropped. Record 5 is missing its minimum age, so this test cannot compare its ages and leaves it unflagged. Record 4 spans 43 million years; no test flags a wide range on its own, and the section on filtering by age range comes back to it.

Record 6 has both coordinates set to zero and fails `zero`, `equal` and `country_mismatch`, since (0, 0) lies in the Gulf of Guinea. Record 7 fails `centroid`: it lies within the default 1 km of the centroid of Germany. Record 8, 3 km from it, passes.

## Filtering records

`keep_unflagged()` keeps the records that no test in `tests` flags. It also checks that every named test was run on the flagged table, and stops with an error if one was not. The later sections reuse this rule as `rule`:

```{r keep}
rule <- tests
kept <- keep_unflagged(f, rule)
kept[, c("occurrence_no", "accepted_name", "min_ma", "max_ma")]
```

Record 5 stays because `age_equal` cannot compare a missing age, so the record carries no flag. A study that needs both ages adds `!is.na(f$min_ma) & !is.na(f$max_ma)` to the selection.

## A wider centroid buffer

OccurrenceCleaner uses default buffers of 1 km around centroids and 100 m around institutions. You can change these with `buffers`.

Fossil localities can be georeferenced coarsely, so a collection placed a few kilometres from a country centroid may have been assigned the country-centroid coordinates. A wider buffer catches such records, at the cost of flagging any correct locality within the same distance. With a 5 km buffer, chosen here to reach record 8, both records near the centroid of Germany are flagged:

```{r buffers}
f5 <- flag_occurrences(fos, tests = tests, columns = cols,
                       buffers = list(centroid = 5000))
data.frame(occurrence_no = fos$occurrence_no,
           centroid_1km = f$flag_centroid, centroid_5km = f5$flag_centroid)
```

These defaults differ from those in CoordinateCleaner. In version 3.0.1, [`clean_fossils()`](https://ropensci.github.io/CoordinateCleaner/reference/clean_fossils.html) applies its defaults `centroids_rad = 0.05` and `inst_rad = 0.001` as radii in metres, so its centroid and institution tests flag only records within 5 cm and 1 mm of a reference point.

## Filtering by age range

Calculate the width of each record's age range as `max_ma - min_ma`. Use `tally_flags()` to count flagged and retained records in each age-range group:

```{r range}
f$age_range <- f$max_ma - f$min_ma
tally_flags(f, bins = list(age_range = c(5, 15)), rules = list(all = rule))
```

Record 4 has an age range of 43 Ma and is the only record above 15 Ma. To retain unflagged records with an age range of at most 15 Ma, apply the limit to the kept records; `which()` also drops records without an age range:

```{r cut}
kept15 <- keep_unflagged(f, rule)
kept15 <- kept15[which(kept15$age_range <= 15), ]
kept15$occurrence_no
```

The 15 Ma limit is illustrative; choose a limit appropriate for your study. You can change it without rerunning the flagging tests.

For comparison, CoordinateCleaner provides this filter through [`cf_range()`](https://ropensci.github.io/CoordinateCleaner/reference/cf_range.html) with `method = "time"`, which flags age ranges wider than `max_range` Ma.

## Outliers within a taxon

Three tests compare each record with others in the same taxon. Here, taxa are the genera in `accepted_name`:

- `age_range_outlier` identifies unusually wide age ranges.

- `age_outlier` identifies unusual ages.

- `spatiotemporal_outlier` identifies unusual combinations of location and age.

The tests use an upper threshold calculated within each taxon: the 75th percentile plus `mltpl` times the interquartile range, with `mltpl = 5` by default. Each test applies this threshold to its own measure of difference. Taxa need at least seven records to be tested.

The two tests that use ages, `age_outlier` and `spatiotemporal_outlier`, sample one age from each record's range. They repeat this five times and flag records identified in a majority of runs. Sampling uses a seed for each taxon, making the results reproducible. You can change the seed through `settings`.

Generate 600 synthetic records with PBDB column names. Locations fall within a region of Wyoming, ages span 46 to 60 Ma, and age ranges are 1 to 4 Ma wide. Genera are assigned independently of location and age.

Then give 20 records a broad age range of 23 to 66 Ma and five records an age range of 99 to 101 Ma. The `injected_case` column labels these changes; unchanged records are labelled `background`.

```{r big}
set.seed(5)
n <- 600
big <- data.frame(
  occurrence_no = 100 + seq_len(n),
  accepted_name = sample(c("Coryphodon", "Hyracotherium", "Mesonyx"), n, TRUE),
  lng = round(runif(n, -110, -106), 3), lat = round(runif(n, 42, 46), 3),
  cc = "US", min_ma = round(runif(n, 46, 56), 1), injected_case = "background")
big$max_ma <- big$min_ma + round(runif(n, 1, 4), 1)
i <- sample(n, 25)
big[i[1:20], c("min_ma", "max_ma", "injected_case")] <- list(23, 66, "broad_range")
big[i[21:25], c("min_ma", "max_ma", "injected_case")] <- list(99, 101, "far_age")
table(big$injected_case)
```

Run the outlier tests alongside the coordinate and age tests to collect their flags in one table:

```{r outliers}
otests <- c("age_range_outlier", "age_outlier", "spatiotemporal_outlier")
g <- flag_occurrences(big, tests = c(tests, otests), columns = cols)
aggregate(as.data.frame(g)[paste0("flag_", otests)],
          by = list(injected_case = g$injected_case), FUN = sum)
```

The range test flags all 20 broad-range records. The age test flags seven broad-range records and all five records near 100 Ma. The spatiotemporal test flags four of the five records near 100 Ma. All three tests leave the background records unflagged.

You can adjust the tests through `settings`. The call below lowers `mltpl` to 2 for the range test and changes the age test's seed to 7:

```{r outlier-settings}
g2 <- flag_occurrences(big, tests = otests, columns = cols,
                       settings = list(age_range_outlier = list(mltpl = 2),
                                       age_outlier = list(seed = 7)))
colSums(as.data.frame(g2)[paste0("flag_", otests)])
```

The range test still flags the same 20 records: their 43 Ma ranges exceed the threshold under both settings. Changing the age test's seed gives 14 flagged records, compared with 12 under the default seed.

## Export to PyRate

PyRate (Silvestro et al. [2014](https://doi.org/10.1093/sysbio/syu006)) estimates speciation and extinction rates from fossil occurrences. Select the records from `g` that are not flagged by any of the tests run and have an age range of at most 15 Ma:

```{r final}
g$age_range <- g$max_ma - g$min_ma
final <- keep_unflagged(g)
final <- final[which(final$age_range <= 15), ]
nrow(final)
```

`export_pyrate()` writes them in PyRate's input format: one array of ages per taxon, an extra age of 0 for extant taxa, and replicates of ages drawn from each record's range.

```{r pyrate}
paths <- export_pyrate(final, file.path(tempdir(), "wyoming"), status = "extinct",
                       columns = cols, replicates = 3, seed = 1)
basename(paths)
```

The `_PyRate.py` file contains the ages and is the input PyRate reads; the `_TaxonList.txt` file lists each taxon with its status. The draws are seeded, so the same call writes the same files; `random = FALSE` uses the midpoints of the ranges instead. CoordinateCleaner's `write_pyrate()` writes the same layout, and the package tests compare the two files.

## Comparison with CoordinateCleaner

`age_equal` and `cf_equal()` give the same flags on the small table:

```{r cc, eval = has_cc}
theirs <- !CoordinateCleaner::cf_equal(fos, value = "flagged", verbose = FALSE)
data.frame(occurrence_no = fos$occurrence_no,
           OccurrenceCleaner = f$flag_age_equal, CoordinateCleaner = theirs)
```

The package's test suite compares the outlier tests with `cf_age()`, `cf_outl()` and `cf_range()`: `age_range_outlier` gives identical flags, and the two tests that draw ages give identical flags when no age has to be drawn.

## Fossils in GBIF data

GBIF stores fossil specimens alongside modern records, marked by `basisOfRecord = "FOSSIL_SPECIMEN"`, usually without PBDB-style age fields. `basis_fossil` flags them:

```{r gbif}
gb <- data.frame(
  gbifID = 1:3,
  basisOfRecord = c("HUMAN_OBSERVATION", "FOSSIL_SPECIMEN", "PRESERVED_SPECIMEN"),
  decimalLongitude = c(16.37, 11.2, 14.5), decimalLatitude = c(48.21, 50.1, 47.1))
fg <- flag_occurrences(gb, tests = c("invalid", "basis_fossil"))
```

A study of present-day distributions keeps the records without the flag, and a palaeontological one keeps only the flagged records:

```{r gbif-select}
keep_unflagged(fg, "basis_fossil")[, c("gbifID", "basisOfRecord")]
fg[has_flag(fg, "basis_fossil"), c("gbifID", "basisOfRecord")]
```

## References

Silvestro, D., Schnitzler, J., Liow, L. H., Antonelli, A. and Salamin, N. (2014). Bayesian estimation of speciation and extinction from incomplete fossil occurrence data. *Systematic Biology* 63:349-367. https://doi.org/10.1093/sysbio/syu006