Introduction

The PREDICTS database (Hudson et al. 2017) is a fully open-access database containing spatial comparisons of biodiversity in different land uses worldwide.

You can get started quickly by downloading the PREDICTS database from my website, but if you intend to use the database in your research, I would ask you please to download the database (freely) from the Natural History Museum data portal, so that appropriate credit is given. The main 2016 database release is here, and an additional 2022 data release here.

For manipulating data, we will use the ‘dplyr’ package and for plotting the ‘ggplot2’ package, which you can install and load as follows:

install.packages("dplyr")
intall.packages("ggplot2")

library(dplyr)
library(ggplot2)

# First we need to download and load the main 2016 release, as well as the 
# supplementary data released in 2022
url("https://www.dropbox.com/scl/fi/0mqoacviqiiurrtplby2m/predicts.rds?rlkey=m5ijge7w1dthkvm5gbp652uel&dl=1") %>% readRDS() %>% droplevels() -> predicts_2016

url("https://www.dropbox.com/scl/fi/3dqv4ptkqkuswxx350xbx/predicts_2022.rds?rlkey=20g6v43l0vjdh3k8e7p533hri&dl=1") %>% readRDS() %>% droplevels() -> predicts_2022

suppressMessages(full_join(predicts_2016,predicts_2022)) -> predicts

To work with the database, you should install a package that contains some basic data-manipulation functions:

remotes::install_github("timnewbold/predicts-demo",subdir="predictsFunctions")

library(predictsFunctions)

Database Structure

The PREDICTS database has a hierarchical structure.

At the highest level, the data in the 2016 release of the PREDICTS database are drawn from 595 published papers (Sources) comparing biodiversity in different land uses.

Each data Source consists of one or more Studies. Sometimes, a published paper will contain biodiversity samples collected using more than one fundamental approach. For example, a paper may contain data for different taxonomic groups that would not appropriately have been sampled using the same methods. In such cases, because data are not directly comparable if collected using different methods, data Sources were split into multiple Studies. Otherwise, each Source contained a single Study.

Each Study consists of data sampled in one or more spatial Blocks. A Study was divided into more than one Block if the sampling showed spatial structuring, to allow us to delimit spatial aggregations of Sites, and thus to account for spatial structuring later when modelling the data.

Each Block contains two or more Sites, which are the specific locations at which biodiversity was sampled.

The complete 2016 version of the PREDICTS database contains data from 817 Studies, which together sample a total of 35,736 Sites.

Within each Site, there are a series of Measurements - of species abundance (3,671,179 records), species occurrence (647,334 records), or species richness (295 records) - for each of the sampled taxa.

Tree map showing the frequency of different biodiversity metric types across the PREDICTS database. Note that counts of metrics were square-root transformed to counter the strong dominance of abundance measures

Key Sampling Information

Measurements in PREDICTS consist of three broad types: species abundance, species occurrence (i.e., presence or absence), or species richness (the number of species sampled). Within these categories, there are lots of different specific metrics. Some measures of diversity are inherently sensitive to sampling effort, which has important implications for processing the data (see below).

Tree map showing the frequency of different biodiversity metric types across the PREDICTS database. Note that counts of metrics were square-root transformed to counter the strong dominance of abundance measures

The different underlying Studies in PREDICTS consist of data sampled using a whole host of different methods, and with widely differing levels of sampling effort. When analysing the PREDICTS data, it is important that diversity Measurements are comparable within Studies, whereas differences among Studies are handled by the hierarchical design of the models (see below). Therefore, it is a requirement of the PREDICTS database that sampling methods are constant within studies. Nevertheless, for 17% of studies, sampling effort differs among the sampled sites. We will see how we deal with this later.

Most of the original biodiversity sampling in the studies contained within the PREDICTS database took place between 2000 and 2013.

Temporal density of mid-point sampling dates for the PREDICTS sites

Sampling at the sites in PREDICTS can span a range of different spatial extents. The maximum linear span of sampling for three-quarters of sites in the PREDICTS database falls between 5 and 490 metres.

Density distribution of maximum linear spans sampled at the PREDICTS sites. The error bar represents the central 75% of values in the distribution

Spatial Information

The vast majority of sampled sites in the PREDICTS database (100%) are associated with geographical coordinates (longitude and latitude).

These geographical coordinates are then used to assign the site within various geographic schemes: biogeographic realms (Ladle and Whittaker 2011; Olson et al. 2001), biomes (Olson et al. 2001), ecoregions (Olson et al. 2001), UN regions and sub-regions (United Nations, 1999), biodiversity hotspots (Myers et al. 2000) and high-biodiversity wilderness areas (Brooks et al. 2006).

The sites in PREDICTS are distributed across all of the world’s 14 terrestrial biomes, 318 ecoregions, and 101 countries.

If you want to reproduce the map below for yourself, you can use the following code, but first we need to install and load some packages for working with map data.

install.packages("rnaturalearth")
intall.packages("tidyterra")

library(rnaturalearth)
library(tidyterra)

suppressMessages(library(sf))

studies <- predicts %>% group_by(SS) %>% 
  summarise(n_spp = n_distinct(Taxon_name_entered),
            Longitude = mean(Longitude),
            Latitude = mean(Latitude))

ne_countries(scale = 50,returnclass = 'sv') -> world

p <- ggplot() + 
  geom_spatvector(data = world,fill="#BF2C2333",col="#00000000") + 
  geom_point(data = studies,
             mapping = aes(x = Longitude,y = Latitude,size=n_spp),
             col = "#1065AB66") +
  scale_x_continuous(name = "") + 
  scale_y_continuous(name = "") + 
  theme_classic()

p

Map of sites with biodiversity samples in the PREDICTS database

Human Pressure Information

The only information on human pressures that is endogenous to the PREDICTS database describes the land-use conditions at the sampled site, with the key variables representing land-use type and intensity of human land use.

In the majority of cases (at least 98.6%), land-use variables in the PREDICTS database were encoded using information given in the original source publications. Occasionally, the data entrant used Google Maps to infer the land-use type.

Land-use type is classified as one of: primary vegetation (apparently pristine habitat, with no record of prior destruction); secondary vegetation (habitat previously destroyed by human actions or extreme natural events, but now recovering towards its natural state); plantation forest (area used for growing woody crops); cropland (area used for growing herbaceous crops); pasture (area used for grazing livestock); urban (human settlements or areas of civic amenity). Secondary vegetation is also divided according to the stage of recovery toward its natural state: young, intermediate or mature. For full definitions, see (Hudson et al. 2017).

Land-use intensity is classified at one of three levels (minimal, light or intense), with the criteria for classification depending on land-use type (see (Hudson et al. 2017) for more details)

Distribution of sites among major land uses (note that all stages of recovery of secondary vegetation have been combined in this figure

Taxonomic Information

The names of the sampled taxa, as given in the original source publications, are entered directly into the database. These names were passed through a taxonomic names resolver to attempt to correct any typographical errors from the original papers. Parsed names were then compared to the Catalogue of Life or manually checked, to try to find the accepted name for the species.

A total of 53,925 taxonomic names are given in the original publications that underpin the PREDICTS database. After taxonomic resolution, a total of 32,489 accepted Latin binomials are identified (many taxa are identified to a higher level than species, or are identified only as morpho-species).

Tree map representing the number of species sampled in different phyla

Data Pre-processing

Before modelling with the PREDICTS data, there are a few basic pre-processing steps that must be carried out. First, for the 0.8% of records that are of an effort-sensitive metric, and where sampling effort differs among sampled sites within the study, we need to correct sampling effort. This is done by rescaling sampling effort within each study, such that the most-sampled site gets a relative effort value of 1, and the sampling effort for all other sites is scaled linearly relative to this maximum value. The CorrectSamplingEffort function in the predictsFunctions package carries out these operations:

# Correct effort-sensitive abundance measures (assumes linear relationship between effort and recorded abundance)
predicts <- predicts %>% CorrectSamplingEffort()

## Correcting 10763 missing sampling effort values
## Rescaling sampling effort
## Correcting 3104852 values for sensitivity to sampling effort

Next, some samples in PREDICTS were entered such that a sampled ‘Site’ consisted of very finely divided samples. For example, samples may have been entered as individual traps. This is not in keeping with the strict definition of a ‘Site’ in the PREDICTS database. To get around this issue, we combine sites with identical coordinates, belonging to the same study and spatial block of sites, sampled on the same dates, sampled with the same methods and recorded using the same diversity metric, and situated within the same land-use type and human land-use intensity. This merging of sites is carried out by the MergeSites function in the predictsFunctions package (warning, this operation can take a few minutes to run, so you might want to go and make a cup of tea!):

predicts <- predicts %>% MergeSites(silent=TRUE)

Basic Models

There are two main ways that you can now analyze the PREDICTS data:

1. Using the raw species-level data, modelling the presence/absence or abundance of individual species; or
2. Using site-level biodiversity summaries, such as species richness, total community abundance or other aggregate biodiversity measures

In this tutorial, we will focus on the second model type, which is the approach used most often in published papers using the PREDICTS database.

We will use the glmmTMB package for modelling.

install.packages("glmmTMB")

library(glmmTMB)

Before we construct some models, we need first to calculate site-level diversity summaries (species richness, total abundance, etc.). This is done using the SiteMetrics function of the predictsFunctions package. We will also rearrange the land-use classification so that sampling is sufficient across different land-use types for creating robust models. To allow the models to run quickly, we will focus on one taxonomic subset of the PREDICTS data: insects (Class ‘Insecta’).

sites <- predicts %>%
  filter(Class=="Insecta") %>%
  droplevels() %>%
  SiteMetrics(extra.cols = c("Longitude","Latitude",
                             "Predominant_land_use",
                             "SSB","SSBS"),
              srEstimators = NULL) %>%
  mutate(LandUse = recode(Predominant_land_use,
                          'Primary vegetation'='Primary',
                          'Mature secondary vegetation'='Secondary_Mature',
                          'Intermediate secondary vegetation'='Secondary_Intermediate',
                          'Young secondary vegetation'='Secondary_Young',
                          'Secondary vegetation (indeterminate age)'=NA_character_,
                          'Plantation forest'='Plantation',
                          'Cropland'='Cropland',
                          'Pasture'='Pasture',
                          'Urban'='Urban',
                          'Cannot decide'=NA_character_)) %>%
  mutate(LandUse = factor(LandUse,levels=c(
    'Primary','Secondary_Mature','Secondary_Intermediate','Secondary_Young',
    'Plantation','Pasture','Cropland','Urban')))

## Computing site metrics for 927934 measurements
## The data contain 244 sources, 344 studies and 8791 sites
## Computing site-level values
## Computing total abundance
## Computing species richness
## Computing Simpson's diversity
## Assembling site-level values

We will start by constructing a very simple model of species richness as a function of land use, fitting the following random intercepts: 1) study identity (SS in the PREDICTS database fields), to control for variation in sampling methods, sampling effort, taxonomic group, and broad study region among individual studies in the PREDICTS database; and 2) Spatial Block nested within Study (SSB in the PREDICTS database fields), to account for the spatial structuring of sampled sites within each study. For now, we will fit a simple Poisson model, but note that this doesn’t conform well to model assumptions because of over-dispersion and zero-inflation in the counts of the number of species in the PREDICTS database. There are more appropriate ways to model the data to handle these issues, but we won’t deal with those yet.

# Exclude NAs from LandUse column, then build model
richMod1 <- sites %>% 
  subset(!is.na(LandUse)) %>%
  glmmTMB(formula=Species_richness~LandUse+(1|SS)+(1|SSB),family=poisson)

To plot the model-estimated effects of land use on insect species richness, we can simply back-transform the coefficients of the model. Since a Poisson model uses a log-transformation on the raw values, we back-transform them by taking the exponential. This tells us what proportion of the number of species in the baseline land use (primary vegetation) are estimated to be retained in the other land-use types, on average. If we multiply these proportions by 100 and subtract 100, this tells us the model-estimated percentage difference in species richness compared to primary vegetation. We can also extract the uncertainty (given as standard error) around these coefficients, from which we can calculate confidence intervals for our modelled estimates. There are more sophisticated ways of plotting model results, but we won’t cover those for now.

# Get the modelled coefficients for the effect of land use on insect species richness
luCoefs <- summary(richMod1)$coefficients$cond[,'Estimate']
# We will replace the first coefficient with zero, since we are going to compare 
# everything to primary vegetation as our baseline
luCoefs[1] <- 0

# Get the uncertainty (standard error) for these land-use coefficients
luCoefSEs <- summary(richMod1)$coefficients$cond[,'Std. Error']
# Again, we will replace the first standard error value with zero, since we are 
# not interested in the uncertainty for our baseline (primary vegetation)
luCoefSEs[1] <- 0

# Now we will calculate the median modelled estimates for the land-use effects,
# taking the exponential to back-transform the raw coefficients from the model
# then multiplying by 100 and substracting 100 to get a percentage difference
y.median <- exp(luCoefs)*100-100
# Now we will calculate the confidence intervals, noting that the 95% confidence
# interval is +/- 1.96 * standard error, and the 67% confidence interval
# +/- 1 * standard error
y.lower.95 <- exp(luCoefs - 1.96*luCoefSEs)*100-100
y.upper.95 <- exp(luCoefs + 1.96*luCoefSEs)*100-100
y.lower.67 <- exp(luCoefs - luCoefSEs)*100-100
y.upper.67 <- exp(luCoefs + luCoefSEs)*100-100
lu <- factor(c("Primary","Secondary - Mature",
        "Secondary - Intermediate","Secondary - Young",
        "Plantation","Pasture","Cropland","Urban"),
        levels=c("Primary","Secondary - Mature",
        "Secondary - Intermediate","Secondary - Young",
        "Plantation","Pasture","Cropland","Urban"))

# Now we can plot the final graph, with median modelled estimates (points) along with
# the confidence intervals (95% as thin lines and 67% as thick lines)
plotLU <- ggplot() + 
  geom_linerange(
    mapping = aes(x = lu,
                  ymin = y.lower.95,
                  ymax = y.upper.95,
                  col = lu)) +
  geom_pointrange(mapping = aes(x = lu,
                                y = y.median,
                                ymin = y.lower.67,
                                ymax = y.upper.67,
                                col = lu),
                  linewidth = 1, size = 0.8) +
  scale_colour_manual(values = c("#66a61e","#147659","#1b9e77","#8ecfbc",
                                 "#7570b3","#d95f02","#e6ab02","#e7298a")) +
  scale_y_continuous(name = "\u0394 Species Richness (%)",
                     limits = c(-30,0)) + 
  scale_x_discrete(name = "Land use") + 
  geom_hline(mapping = aes(yintercept = 0),
             alpha = 0.3, linetype = 'dashed') +
  theme_classic() + 
  theme(legend.position = "none") + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

plotLU

Model estimated effects of land use on insect species richness. Points show median predicted species richness as a % difference compared to primary vegetation, thick lines show 67% credible intervals, and thin lines 95% credible intervals