Worked example 3: What does a near-random dataset look like?

if (!"occAssess" %in% installed.packages()) devtools::install_github("https://github.com/robboyd/occAssess")
library(occAssess)

Introduction

This vignette provides a worked example for the functionality of occAssess. For a full tutorial see the main vignette supplied in Boyd et al. (2021).

Occurrence data

In this worked example, we simulated data that approximates random distributions in space and time for fourty species. The species were simulated with varying prevalence (number of records varies by a factor of 5 with a mean of ~1000). The data were generated at random locations in the UK, and uniformly over the period 2001 to 2010. Each simulated data point was randomly assigned to one of two surveys which are specified in the “identifier” field. The data can be accessed within occAssess as follows:


data(random40Species)

spDat <- random40Species

str(spDat)
#> 'data.frame':    40000 obs. of  25 variables:
#>  $ species           : Factor w/ 40 levels "species1","species10",..: 14 19 8 28 5 33 12 27 36 12 ...
#>  $ x                 : num  213500 289500 369500 223500 374500 ...
#>  $ y                 : num  797500 745500 537500 591500 520500 ...
#>  $ year              : int  2009 2004 2006 2010 2009 2004 2004 2004 2010 2001 ...
#>  $ spatialUncertainty: num  9802 3506 15310 9810 7189 ...
#>  $ identifier        : Factor w/ 2 levels "survey1","survey2": 2 2 1 1 1 1 1 1 1 1 ...
#>  $ clim1             : num  7.01 6.65 6.76 8.05 7.3 ...
#>  $ clim2             : num  12.2 12.2 12.4 13 12.9 ...
#>  $ clim3             : num  2.43 1.82 1.82 3.63 2.32 ...
#>  $ clim4             : num  2402 1300 1464 1364 1368 ...
#>  $ clim5             : num  296 167 164 152 149 ...
#>  $ clim6             : num  105.6 67.9 84.6 70.7 80.4 ...
#>  $ clim7             : num  36.1 29 24.6 26.2 23.6 ...
#>  $ clim8             : num  862 440 483 446 445 ...
#>  $ clim9             : num  341 219 259 229 245 ...
#>  $ clim10            : num  399 240 292 280 282 ...
#>  $ clim11            : num  780 412 440 392 400 ...
#>  $ clim12            : num  6.39 6.79 6.87 6.2 6.92 ...
#>  $ clim13            : num  38 37.5 37.1 37.5 37.1 ...
#>  $ clim14            : num  406 430 435 392 437 ...
#>  $ clim15            : num  16.5 16.9 17.5 17.2 18.1 ...
#>  $ clim16            : num  -0.322 -1.163 -1.018 0.666 -0.573 ...
#>  $ clim17            : num  16.8 18.1 18.5 16.5 18.6 ...
#>  $ clim18            : num  3.11 2.55 2.6 6.3 3.12 ...
#>  $ clim19            : num  8.48 8.2 8.21 9.29 8.75 ...

Periods

In this example, we will specify five periods over 2001 to 2010


periods <- list(2001:2002, 2003:2004, 2005:2006, 2007:2008, 2009:2010)

Functions

All of the functions in occAssess require two common arguments: dat and periods (outlined above). I will run through each function in the following, indicating where additional arguments are required. Generally, the functions in occAssess return a list with two elements: one being a ggplot2 object, with a separate panel for each level of identifier; and a second with the data underpinning the plot.

assessRecordNumber

The first function I will introduce is the simplest: assessRecordNumber. This function simply plots out the number of records per year in your dataset.

This function enables researchers to quickly establish how the number of records has changed over time.

assessSpeciesNumber

In addition to the number of records, you may wish to know how the number of species (taxonomic coverage) in your dataset changes over time. For this you can use the function assessSpeciesNumber:

assessSpeciesID

It has been speculated that apparent changes in taxonomic coverage could, in fact, reflect a change in taxonomic expertise over time. For example, if fewer individuals have the skill to identify certain species, then it may not appear in your dataset in the later periods. The function assessSpeciesID treats the proportion of species identified to species level as a proxy for taxonomic expertise:

The argument “type” can take the values proportion (proportion of records identified to species level) or count (number of records identified to species level).

assessRarityBias

A number of studies have defined taxonomic bias in a dataset as the degree of proportionality between species’ range sizes (usually proxied by the number of grid cells on which it has been recorded) and the total number of records. One can regress the number of records on range size, and the residuals give an index of how over-or undersampled a species is given its prevalence. The function assessRarityBias conducts these analyses for each time period, and uses the r2 value from the linear regressions as an index proportionality between range sizes and number of records. Higher values indicate that species’ are sampled in proportion to their range sizes whereas lower values indicate that some species are over- or undersampled.

Note the warning message which tells us that there are low numbers of species in some periods (not surprising as the data only contain five species). This represents a small sample size for the regression of range size on number of records so the results should be viewed with caution.

assessSpatialCov

The function assessSpatialCov grids your data at a specified spatial resolution then maps it in geographic space. In this example, I provide a shapefile with the boundaries of the UK. If I was working on the WGS84 coordinate reference system (here I am using OSGB 36) this would not be necessary; I could instead use the countries argument and simply specify “UK”.

As you can see there are three new arguments to be specified. res is the spatial resolution at which you would like to map the data (units depend on you coordinate reference system, e.g. m if easting and northing, and decimal degress in lon/ lat); logCount indicates whether or not you would like to log10 transform the counts for visual purposes; and countries defines the countries covered by your data. Countries must be specified in order to plot their boundaries.

assessSpatialUncertainty

Point occurrence data often comes with assocatiated spatial uncertainty (i.e. how uncertain the coordinates are in x and y dimensions). For example, GBIF data comes with a field called coorinateUncertaintyInMeters. The function assessSptialUncertainty can be used to visualize spatial uncertainty in your dataset as a kernel density plot:

assessSpatialBias

Even if your data has good spatial coverage, it may be biased; that is to say, it may deviate from a random distribution in space. The function assessSpatialBias provides an index of how far your data deviates from a random distribution. To do this is simulates an equal number of points to your data randomly across your study region. Then, for each time period, it calculates the average nearest neighbour distance across your data points and divides it by the average nearest neighbour distance from the random sample. If the index is lower than one then your data is more clustered than the random sample, and if it is above one it is more dispersed. To delineate your study area, you must provide a mask layer. The mask is a raster object which is has numeric values within your study area, and is NA outside of your study area. Here, I’ll use some species distribution model outputs for the UK as a mask layer:

The argument nSamps indicates how many random distributions should be drawn, and the argument degrade = TRUE indicates that any duplicated coordinates within a time period and for a given level of identifier are removed. The shaded regions on the plot indicate the 5th and 95th percentiles of the nearest neighbour index calculated over nSamps random samples.

assessEnvBias

Spatial bias in your dataset does not necessarily tell you anything about environmental bias. The function assessEnvBias assess the degree to which your data are biased across time periods in environmental space. To do this we first need to get some climate data. I will use the standard suite of 19 bioclimatic variables from worldclim. It is possible to get this data through R using the raster package, but here I will use my local version for speed:

assessEnvBias conducts a principal component analysis on your environmental data, then maps your occurrence data in environmental space:

It is also possible to modify the appearance of envBias$plot using additional arguments that can be passed to ggfortify::autoplot. For example, you can include elipses, use different principal components, include vaiable vectors, etc.