Spatial subsampling tutorial

Contents

• Why to subsample spatially

• Prepare Paleobiology Database records

• Subsampling examples

  • Circular subsampling

  • PBDB collections, references, and duplicate entries

  • Nearest-neighbour subsampling

  • Latitudinal subsampling

• Analysis on subsampled data

• References

This document introduces rationale for spatial subsampling and covers common use cases for divvy functions, including example R code. Sections below discuss ways to format data, implement three different subsampling routines, and calculate common biodiversity metrics.

For a more comprehensive review of spatial subsampling theory, applications in paleobiology and ecology, and implementation in the divvy package, refer to Antell, Benson, and Saupe (accepted), ‘Spatial standardization of taxon occurrence data—a call to action’.

Why to subsample spatially

Biodiversity measures, whether quantified from fossil or extant data, reflect the outcome of three types of processes: 1. the biological processes the measurements aim to capture, 2. stochastic processes that generate random error, and 3. observation processes, such as the distribution of collecting effort intensity. Compared to neontological data, observation densities of palaeontological data reflect the additional contributions of taphonomy (decay and preservation processes) and sedimentation (burial and lithification processes). The term ‘sampling’ is used herein to refer to the full suite of influences on how many of the fossil taxa that lived in a place are recorded as fossil occurrences there. Given the large number of processes involved in sampling fossils, it is unsurprising that for many palaeo-occurrence datasets, the combined effect size of unexplained error and variable observation probabilities on biodiversity parameter estimates is greater than that of true biological processes from the past.

The relative strength of sampling structure compared to biological structure is a well-known issue in quantitative palaeobiology and therefore has received considerable attention. The most mainstream method of standardising sampling in biodiversity data (fossil or extant) is rarefaction, a genre of statistical method that resamples taxon occurrences to a given threshold of sampling completeness. ‘Classical’ rarefaction uses cumulative specimen count as a proxy of sampling completeness, while ‘coverage-based’ rarefaction estimates the cumulative coverage of a taxon frequency distribution curve as a proxy (Alroy 2010; Chao and Jost 2012). Examples of both types of rarefaction abound in recent palaeobiology literature. Rarefaction largely succeeds in producing diversity estimates corrected for sampling differences between sites, i.e. diversity estimates at discrete sites (alpha diversity) that can be compared fairly. However, challenges arise in comparing diversity between times or regions that contain multiple sites.

Taxonomic composition turns over from one geographic location to another. The amount of this turnover (beta diversity) generally increases with distance, a relationship known as the species–area effect. In practical terms, the greater the spatial dispersion or area of sampling, the more taxa will tend to be recovered. Sampling processes control both the dispersion and area of localities with fossil data, which makes it essential for palaeobiologists to account for spatial coverage of sampling when conducting biodiversity analyses, regardless of the use of rarefaction (Benson et al. 2021). Without standardising sampling in a spatial context, any estimates of biodiversity parameters will be an indiscernible mix of true historic values conflated with the amount of observation of those values.

divvy implements several spatial subsampling methods that control both the geographic dispersion and area of taxon occurrences, so analysts can make fairer estimates of biodiversity parameters between regions or through time. Each method is iterative, drawing many, comparable spatial replicates.

Prepare Paleobiology Database records

We begin by loading the divvy package itself, and note its geospatial package dependencies units, sf, terra, and vegan, and the iNEXT package for taxonomic richness rarefaction.

library(divvy) 
library(units)
library(sf)
library(terra)
library(vegan)
library(iNEXT)

divvy includes the attached dataset bivalves, a download of fossil occurrences from the Paleobiology Database (PBDB) that has been subset to a few relevant columns and minimally cleaned. bivalves contains ca. 8,000 (palaeo)coordinates and identifications for ca. 500 marine bivalve genera from the Pliocene. For more information about the dataset, query ?bivalves.

The attached dataset is provided only as an example for working with PBDB-structured occurrence records; formal analyses would be wise to vet downloaded data more rigorously, including revising taxonomy, time bin assignments, and environmental classifications. The fossilbrush package provides a palette of tools to help clean PBDB data.

data(bivalves)
head(bivalves)

#>         genus paleolng paleolat collection_no reference_no
#> 1        Lima   176.85   -42.69         34849         9315
#> 2 Nemocardium   176.85   -42.69         34849        11706
#> 3      Ostrea   176.85   -42.69         34849        11706
#> 4     Chlamys   176.85   -42.69         34849         9315
#> 5  Pycnodonte  -114.37    33.76         34857         9338
#> 6      Euvola  -114.37    33.76         34857         9338
#>              environment max_ma min_ma          accepted_name
#> 1             basin reef  5.333    3.6                   Lima
#> 2             basin reef  5.333    3.6 Nemocardium (Pratulum)
#> 3             basin reef  5.333    3.6       Ostrea chilensis
#> 4             basin reef  5.333    3.6                Chlamys
#> 5 marginal marine indet.  5.333    3.6   Pycnodonte heermanni
#> 6 marginal marine indet.  5.333    3.6           Euvola keepi

nrow(bivalves)

#> [1] 8095

length(unique(bivalves$genus))

#> [1] 550

The latitude-longitude coordinates in the PBDB come from many different studies, which means the precision and accuracy vary across records. Fossils collected from the same locality may be reported with different coordinates due purely to different GPS equipment, mapping, or decimal rounding, for example. The spatial subsampling procedures below put great weight on the number of unique localities in a region, so it is important to avoid inflating the site number artificially.

One standard way to smooth over slight discrepancies in occurrence positions is to convert point coordinates (‘vector’ spatial data) to grid cells (‘raster’ spatial data)¹. Effectively, one lays a web over a study region and records the position of polygons (‘grid cells’) in which points fall, rather than the original xy point coordinates. An additional advantage of raster over vector data is the tremendous increase in efficiency of spatial computations.

Palaeobiologists often convert abundance counts to binary presence-absence data for taxon occurrences in grid cells. This practice is especially common for PBDB data because abundance information is usually non-standard, if it is recorded at all. (Read below for further notes about duplicate taxon occurrences.)

Many PBDB analyses are also global in extent, which makes the choice of coordinate reference system for the raster grid important. The areas and shapes of grid cells at the poles vs. the equator can differ widely with certain reference systems or map projections. The code below converts latitude-longitude coordinates to the Equal Earth reference system, an equal-area projection for world maps. Another option of friendly raster grid system is the is the icosa package, developed by palaeobiologist Ádám Kocsis, which creates tessellations of pentagons and hexagons. The polygons are approximately equal in area and shape across the globe. The spatial resolution of the grid is controlled by arguments in the hexagrid function.

# initialise Equal Earth projected coordinates
rWorld <- rast()
prj <- 'EPSG:8857'
rPrj <- project(rWorld, prj, res = 200000) # 200,000m is approximately 2 degrees
values(rPrj) <- 1:ncell(rPrj)

# coordinate column names for the current and target coordinate reference system
xyCartes <- c('paleolng','paleolat')
xyCell   <- c('cellX','cellY')

# retrieve coordinates of raster cell centroids
llOccs <- vect(bivalves, geom = xyCartes, crs = 'epsg:4326')
prjOccs <- project(llOccs, prj)
bivalves$cell <- cells(rPrj, prjOccs)[,'cell']
bivalves[, xyCell] <- xyFromCell(rPrj, bivalves$cell)

We can inspect the spatial distribution of data by converting occurrence points to spatial features (a class that contains coordinate system information) and plotting them against a basic world map, supplied here from the rnaturalearth and rnaturalearthdata packages.

occUniq <- uniqify(bivalves, xyCell)
ptsUniq <- st_as_sf(occUniq, coords = xyCell, crs = prj)

library(rnaturalearth)
library(rnaturalearthdata)
library(ggplot2) # for plot visualisation
world <- ne_countries(scale = "medium", returnclass = "sf")
worldP <- ggplot(data = world) +
  theme_bw() +
  geom_sf() +
  geom_sf(data = ptsUniq, shape = 17, color = 'blue')

plot(worldP)

Subsampling examples

divvy offers three approaches to spatially subsample data:

• cookies: Imposes a radial constraint on the spatial bounds of a subsample and standardises area by rarefying the number of localities

• clustr: Aggregates sites that are nearest neighbours (connecting them with a minimum spanning tree) to impose a maximum diameter on the spatial bounds of a subsample, and optionally rarefies localities

• bandit: Rarefies the number of localities within bands of equal latitude

set.seed(1)
circLocs <- cookies(dat = bivalves,  
                    xy = xyCell,
                    iter = 500, 
                    nSite = 12, 
                    r = 1500, # radial distance in km
                    weight = TRUE, # probabilistically aggregate subsampling sites
                    crs = prj, # Equal Earth projection
                    output = 'locs')
length(circLocs)

#> [1] 500

circLocs[[1]]

#>         cellX   cellY
#> 1019 -6543959 4739440
#> 56   -7343959 3539440
#> 7586 -7143959 3939440
#> 1278 -7343959 3939440
#> 1049 -6543959 4539440
#> 2921 -7143959 3539440
#> 1054 -6743959 4339440
#> 1076 -6743959 4539440
#> 122  -7743959 3939440
#> 1251 -7543959 3939440
#> 2065 -7543959 3739440
#> 911  -6543959 4339440

# over-plot the subsample locations
smplPts <- st_as_sf(circLocs[[1]], coords = xyCell, crs = prj)
worldP +
  geom_sf(data = smplPts, shape = 17, color = 'red')

cntr <- smplPts[1,]
# distances inferred to be meters based on lat-long coord system
r <- 1500
buf <- st_buffer(cntr, dist = r*1000)

# over-plot subsample boundary region and included sites
worldP +
  geom_sf(data = smplPts, shape = 17, color = 'red') +
  geom_sf(data = buf, fill = NA, linewidth = 1, color = 'red')

inBuf <- st_intersects(ptsUniq, buf, sparse = FALSE) 
sum(inBuf)

#> [1] 14

set.seed(7)
# same parameter values as above except for 'output'
circOccs <- cookies(dat = bivalves, 
                    xy = xyCell, 
                    iter = 500, 
                    nSite = 12, 
                    r = 1500, 
                    weight = TRUE, 
                    crs = prj,
                    output = 'full')
head( circOccs[[8]] )

#>           genus paleolng paleolat collection_no reference_no    environment
#> 83      Anadara   -69.82    10.25         41552        11132  marine indet.
#> 409  Iliochione   -79.18     0.54         51887        13814       offshore
#> 410 Undulostrea   -79.18     0.54         51887        13814       offshore
#> 411  Argopecten   -79.18     0.54         51887        13814       offshore
#> 412  Argopecten   -79.18     0.54         51881        13814 coastal indet.
#> 413        Arca   -79.18     0.54         51881        13814 coastal indet.
#>     max_ma min_ma                accepted_name cell    cellX     cellY
#> 83   5.333  3.600 Anadara (Larkinia) notabilis 6074 -6543959 1339439.8
#> 409  5.333  2.588      Iliochione callistoides 7101 -7543959  139439.8
#> 410  5.333  2.588          Undulostrea megodon 7101 -7543959  139439.8
#> 411  5.333  2.588       Argopecten ventricosus 7101 -7543959  139439.8
#> 412  5.333  2.588       Argopecten ventricosus 7101 -7543959  139439.8
#> 413  5.333  2.588                         Arca 7101 -7543959  139439.8

bivUniq <- uniqify(bivalves, taxVar = 'genus', xyCell)
nrow(bivUniq)

#> [1] 3061

set.seed(2)
nnLocs <- clustr(dat = bivalves, 
                 xy = xyCell, 
                 iter = 500, 
                 distMax = 3000, # diameter = 2x the circular radius set above
                 nSite = 12,
                 crs = prj
                 )
nnLocs[[1]]

#>         cellX      cellY
#> 81   -6343959 1339439.78
#> 1891 -7343959  -60560.22
#> 5884 -7343959  539439.78
#> 6477 -7743959 -460560.22
#> 1491 -5743959 1339439.78
#> 2415 -7943959 1139439.78
#> 409  -7543959  139439.78
#> 5916 -6543959 1539439.78
#> 1936 -6743959 2339439.78
#> 3887 -7743959 1339439.78
#> 74   -5943959 1339439.78
#> 5751 -7143959 1339439.78

set.seed(2)
nnAllSites <- clustr(dat = bivalves, 
                     xy = xyCell, 
                     iter = 500,
                     distMax = 3000, # diameter = 2x the circular radius set above
                     nMin = 3,
                     crs = prj
                     )
nrow( nnAllSites[[1]] )

#> [1] 17

bandLocs <- bandit(dat = bivalves,
                   xy = xyCell,
                   iter = 100, nSite = 12, 
                   bin = 20, # interval width in degrees
                   crs = prj
                   # ,absLat = TRUE # optional
                   )
nrow(bandLocs[[1]]) # number of sites in one subsampled band

#> [1] 12

length(bandLocs) # number of subsamples, tallied across all bands

#> [1] 500

unique(names(bandLocs)) # latitudinal degrees of subsampled intervals

#> [1] "[-50,-30)" "[-10,10)"  "[10,30)"   "[30,50)"   "[50,70)"

Analysis on subsampled data

The sdSumry() function returns a summary of the spatial characteristics of a dataset/subsample: number of unique locations, centroid coordinates, latitudinal range (degrees), great circle distance (km), and summed minimum spanning tree length (km) for occurrences. sdSumry() also tallies taxa in a sample and performs coverage-based rarefaction (if quotaQ supplied) and classical rarefaction (if quotaN supplied). Rarefied estimates are returned along with their associated 95% confidence interval (estimated by iNEXT).

When coverage-based rarefaction is applied, the specified coverage level (quotaQ) should be appropriate for the estimated coverage of the original sample. Any requested quota greater than empirically available will require diversity extrapolation. sdSumry() returns the estimated sample coverage whenever coverage-based rarefaction is applied. Also note, homogeneous evenness is an underlying assumption for fair comparisons of coverage-based rarefaction diversity estimates. To help check this assumption, sdSumry() also returns Pielou’s J evenness metric whenever coverage-based rarefaction is applied. Evenness metrics are an area of ongoing theoretical debate and methodological development in ecology, and the use of Pielou’s J in divvy is a reflection of this metrics’ widespread use rather than a singular endorsement. If J estimates vary widely between time intervals or other categories of comparison, it is worth investigating evenness differences in more nuanced detail, for instance by calculating a variety of metrics.

Compare the summary data from the original dataset vs. a single spatial subsample. First consider the summary of spatial metrics. The geographic dispersion of bivalves occurrences is enormous—this in itself indicates spatial standardisation is necessary to derive meaningful biological metrics. When we compare diversity estimates from the global vs. spatially-subsampled dataset, richness estimates are much higher for the global dataset, regardless of rarefaction method. This result demonstrates the species–area effect described in the introduction: even when rarefaction is applied to the same quota or coverage level, diversity estimates are larger for datasets with greater spatial coverage, because rarefaction fails to control spatial turnover in diversity.

unsamp <- sdSumry(dat = bivalves, 
                  taxVar = 'genus',
                  collections = 'collection_no',
                  xy = xyCell,
                  quotaQ = 0.8, quotaN = 100,
                  omitDom = TRUE,
                  crs = prj)
unsamp

#>   nColl nOcc nLoc centroidX centroidY latRange greatCircDist meanPairDist
#> 1   863 3061  157 -434404.9   1647720 137.4159      28917.12     11375.07
#>   minSpanTree     SCOR nTax  evenness coverage  SQSdiv SQSlow95 SQSupr95
#> 1    93349.91 20.64401  550 0.9061343   0.9414 323.015 310.6638 335.3662
#>      CRdiv  CRlow95  CRupr95
#> 1 470.9424 454.3334 487.5515

samp1 <- sdSumry(dat = circOccs[[1]], 
                 taxVar = 'genus',
                 collections = 'collection_no',
                 xy = xyCell,
                 quotaQ = 0.8, quotaN = 100, 
                 omitDom = TRUE,
                 crs = prj)
samp1

#>   nColl nOcc nLoc centroidX centroidY latRange greatCircDist meanPairDist
#> 1   134  442   12  -7027292   4139440 10.25901      1612.452     779.4877
#>   minSpanTree     SCOR nTax  evenness coverage  SQSdiv SQSlow95 SQSupr95
#> 1    2648.528 46.01033  164 0.9572183   0.8704 136.402 123.5963 149.2078
#>      CRdiv CRlow95  CRupr95
#> 1 219.8183 190.056 249.5806

Iterative spatial subsampling addresses this issue in that it can control for beta diversity in global occurrence datasets. Even though any individual subsample includes only a fraction of the total occurrences, all data can contribute to analyses through repeated subsampling. Ecological variables calculated on any given subsample are directly comparable to those from other subsamples on the same dataset or a comparison dataset (e.g. a different time step or habitat type). The code chunk below demonstrates how to calculate summary statistics over all subsamples of Pliocene bivalve occurrences. The median taxon count in a 1500km-radius subsample is 162, with an interquartile range of 126–183 taxa. This range can be interpreted loosely as primary regional variation in estimated palaeodiversity.

# warning - it's slow to rarefy diversity for hundreds of subsamples!
# this code chunk skips it for quick demonstration purposes
sampsMeta <- sdSumry(dat = circOccs, 
                     taxVar = 'genus',
                     collections = 'collection_no',
                     xy = xyCell,
                     crs = prj
                     )
quantile(sampsMeta$nTax, c(0.25, 0.5, 0.75))

#> 25% 50% 75% 
#> 126 162 183