RStudio project

Open the RStudio project that we created in the first session. I recommend to use this RStudio project for the entire course and within the RStudio project create separate R scripts for each session.

  • Create a new empty R script by going to the tab “File”, select “New File” and then “R script”
  • In the new R script, type # Session b5: Pseudo-absence and background data and save the file in your folder “scripts” within your project folder, e.g. as “b5_PseudoAbsences.R”

1 Introduction

In the previous sessions, we have worked with very convenient presence/absence data to train our species distribution models (SDMs). However, as you have seen when downloading your own GBIF data, we often have only presence records available. Absence data are also inherently difficult to get because it requires a very high sampling effort to classify a species as absent from a location. For plants, we would need complete inventories of a larger region, e.g. several 100 m plots placed within larger 1 km sample squares. For birds, we would need several visits to a specific region within the breeding season. In both cases, we could still miss some species, meaning that we do not record the species although present (resulting in a false negative).

But what should we do if absence data are not available? Most SDM algorithms (except the profile methods, see session 6) need some kind of absence or background data to contrast the presence data to. There are different approaches for creating background data or pseudo-absence data (Barbet-Massin et al. 2012; Kramer-Schadt et al. 2013; Phillips et al. 2009), although there is still room for further developments in this field and more clear-cut recommendations for users would be certainly useful. Nevertheless, I hope this tutorial will provide some examples of how to deal with presence-only data. For advice on how many background/pseudo-absence points you need, please read (Barbet-Massin et al. 2012).

In this session, we look at two major ways of creating background/pseudo-absence data:

  • random selection of points within study area (including or excluding the presence locations) (Barbet-Massin et al. 2012)
  • random selection of points outside of study area (Barbet-Massin et al. 2012)

We will not cover approaches for dealing with sampling bias:

  • accounting for spatial sampling bias using target-group selection (Phillips et al. 2009)
  • accounting for spatial sampling bias using inverse distance weighting (Kramer-Schadt et al. 2013)

2 Species presence data

I created a virtual species data set with presence points for a dummy species called Populus imagines and sister species Populus spp. The spatial resolution of the data is 5 minutes. You can download the data here or from the moodle page.

Load and plot the data:


# Our study region and species data
region <- terra::rast('data/Prac5_Europe5min.grd') # mask of study region
sp <- read.table('data/Prac5_presences.txt', header=T) # species presences

# Plot the map and data

3 Background/Pseudo-absence data selection

We use a lot of methods from the predicts and terra tutorials, which are worth looking into (see link here).

3.1 Random selection of points within study area but excluding the presence location

The terra package has a function spatSample() to sample random points (background data).

# Randomly select background points from study region
# The argument na.rm=T ensures that we only sample points within the masked regions (not in the ocean)
# The argument as.points=T indicates that a SpatVector of coordinates should be returned
bg_rand <- terra::spatSample(region, 500, "random", na.rm=T, as.points=TRUE)

# Plot the map and data

By default the function spatSample() will sample from the entire study area independent of the presence points - corresponding to the idea of background data. In order not to sample pseudo-absence at presence locations, we have to tweak our regional mask and assign NA to the presence location

# Make a new regional mask that contains NAs in presence locations:
sp_cells <- terra::extract(region, sp[sp$sp=='Populus_imagines',1:2], cells=T)$cell
region_exclp <- region
values(region_exclp)[sp_cells] <- NA

# Randomly select background data but excluding presence locations
bg_rand_exclp <- terra::spatSample(region_exclp, 500, "random", na.rm=T, as.points=TRUE)

# Plot the map and data

Also, we can define an extent from where random points should be drawn.

# Define extent object:
e <- ext(8,24,46,57)

# Randomly select background data within a restricted extent and excluding presence locations:
bg_rand_exclp_e <- terra::spatSample(region_exclp, 500, "random", na.rm=T, as.points=TRUE, ext=e)

# Plot the map and data
lines(e, col='red')

Last, we could also restrict the random samples to within a certain buffer distance. For this, we first create a SpatVector, then place a buffer around these and sample from within the buffer. The buffer can be interpreted as the area that was accessible to the species in the long-term past. Thus, the buffer helps accounting for biogeographic history by constraining background/pseudo-absence points to those geographic area that could have been reached by the species given the movement capacity but excludes other areas. Of course, restricting the absences to a certain geographic rectangle can achieve a similar tasks but might be less accurate for complex geographies and large areas. For example, is Iceland accessible to European species or not?

# Create SpatVector object of known occurrences:
pop_imag <- terra::vect( as.matrix( sp[sp$sp=='Populus_imagines',1:2]) , crs=crs(region))

# Then, place a buffer of 200 km radius around our presence points
v_buf <- terra::buffer(pop_imag, width=200000)

# Set all raster cells outside the buffer to NA
region_buf <- terra::mask(region, v_buf)

# Randomly select background data within the buffer
bg_rand_buf <- terra::spatSample(region_buf, 500, "random", na.rm=T, as.points=TRUE)

# Plot the map and data
plot(region_buf, legend=F, add=T)

3.2 Random selection of points outside of study area

Barbet-Massin et al. (2012) also suggested a method to sample pseudo-absences only beyond a minimum distance from the presence points. This is more of a macroecological approach, suitable for characterising the climatic limits of species. We can also use the buffering approach from above to achieve this.

# Place a buffer of 200 km radius around our presence points
v_buf <- terra::buffer(pop_imag, width=200000)

# Set all raster cells inside the buffer to NA using the argument inverse=T
region_outbuf <- terra::mask(region, v_buf, inverse=T)
region_buf <- terra::mask(region, v_buf)  # we use this only for illustrative purpose below

# Randomly select background data outside the buffer
bg_rand_outbuf <- terra::spatSample(region_outbuf, 500, "random", na.rm=T, as.points=TRUE)

# Plot the map and data
plot(region_buf,legend=F, add=T)

4 Workflow for joining presence and background/pseudo-absence data

In session b2, we have already learned how to join species data and environmental data (at 10 min resolution) using the Alpine Shrew as example species. Here, we will revisit this example, sample random background data and join these with environmental data to come up with a dataset containing presences and background data as well as the climatic predictors.

# Load the species presence data (here, the data set from session 1):

##          key             scientificName decimalLatitude decimalLongitude
## 1 4416901401 Sorex alpinus Schinz, 1837        47.03705         9.045486
## 2 4156615500 Sorex alpinus Schinz, 1837        48.02569         8.125231
## 3 4156625877 Sorex alpinus Schinz, 1837        48.02580         8.124819
## 4 4416848466 Sorex alpinus Schinz, 1837        46.96636         9.000610
## 5 4129978783 Sorex alpinus Schinz, 1837        47.05048        11.417289
## 6 4133676752 Sorex alpinus Schinz, 1837        47.44182        13.660217
# Plot the species presences

maps::map('world',xlim=c(-12,45), ylim=c(35,73))
points(gbif_shrew_cleaned$decimalLongitude, gbif_shrew_cleaned$decimalLatitude, col='red',  pch=19)

You have already downloaded the climate data at 10-min spatial resolution in session b2 and can simply read it back in. We crop it to European extent.


# Download global bioclimatic data from worldclim (you may have to set argument 'download=T' for first download, if 'download=F' it will attempt to read from file):
clim <- geodata::worldclim_global(var = 'bio', res = 10, download = F, path = 'data')

# Crop to Europe
clim <- terra::crop(clim, c(-12,45,35,73))

What is important to consider is that you kind of arbitrarily chose a scale of analysis by chosing climate data (or other environmental data) at a certain spatial resolution. Here, we chose a spatial resolution of 10 minutes while the species data may actually be at a finer resolution. So, we first make sure that our species data are fit to the spatial resolution, meaning we remove any duplicates within 10 minute cells. We do this by joining the species data with the environmental data (basically, repeating what we had already done at the end of session 2). Then, we can remove the duplicate cells.

# We have already extracted environmental data and raster cellnumbers for the species data
# Remember to remove any rows with duplicate cellnumbers if necessary:
## logical(0)
# From now on, we just need the coordinates:
gbif_shrew_coords <- gbif_shrew_cleaned[,c('decimalLongitude','decimalLatitude')]

We place a buffer of 200 km around the shrew records and sample background points randomly from within the buffer but excluding presence locations (remember that also other pseudo-absence/background data strategies are possible).

# Make SpatVector:
presences <- terra::vect( as.matrix(gbif_shrew_coords), crs=crs(clim))

# Then, place a buffer of 200 km radius around our presence points
v_buf <- terra::buffer(presences, width=200000)

# Create a background mask with target resolution and extent from climate layers
# Set all raster cells outside the buffer to NA.
bg <- clim[[1]]
values(bg)[!] <- 1
region_buf <- terra::mask(bg, v_buf)
plot(bg, col='grey90', legend=F)
plot(region_buf, add=T, col='grey60', legend=F)

# Exclude presence locations:
sp_cells <- terra::extract(region_buf, presences, cells=T)$cell
region_buf_exclp <- region_buf
values(region_buf_exclp)[sp_cells] <- NA

# Randomly select background data within the buffer, excluding presence locations. We sample 10 times as many background data as we have presences. To ensure that we find enough samples, we use the argument exhaustive=T
bg_rand_buf <- terra::spatSample(region_buf_exclp, length(presences)*10, "random", na.rm=T, as.points=TRUE, exhaustive=T)
## Warning: [spatSample] fewer samples than requested are available
## Warning: [spatSample] fewer cells returned than requested
points(bg_rand_buf, pch=19, cex=0.2)
points(presences, pch=19, cex=0.5, col='red')

Next, we need to join the presences and background data, and extract the environmental data.

# First, we prepare the presences data to contain a column indicating 1 for presence.
sp_env <- data.frame(gbif_shrew_coords, occ=1)

# Second, we make sure the background data have the same columns, and indicate 0 for absence.
bg_rand_buf_df <- data.frame(terra::geom(bg_rand_buf)[,c('x','y')])
##        x               y        
##  Min.   : 2.75   Min.   :40.92  
##  1st Qu.: 7.75   1st Qu.:45.08  
##  Median :11.25   Median :47.08  
##  Mean   :10.92   Mean   :47.06  
##  3rd Qu.:14.08   3rd Qu.:49.08  
##  Max.   :18.58   Max.   :52.58
names(bg_rand_buf_df) <- c('decimalLongitude','decimalLatitude')
bg_rand_buf_df$occ <- 0
##  decimalLongitude decimalLatitude      occ   
##  Min.   : 2.75    Min.   :40.92   Min.   :0  
##  1st Qu.: 7.75    1st Qu.:45.08   1st Qu.:0  
##  Median :11.25    Median :47.08   Median :0  
##  Mean   :10.92    Mean   :47.06   Mean   :0  
##  3rd Qu.:14.08    3rd Qu.:49.08   3rd Qu.:0  
##  Max.   :18.58    Max.   :52.58   Max.   :0
# Third, we bind these two data sets
sp_env <- rbind(sp_env, bg_rand_buf_df)
##  decimalLongitude decimalLatitude      occ         
##  Min.   : 2.750   Min.   :40.92   Min.   :0.00000  
##  1st Qu.: 7.583   1st Qu.:45.25   1st Qu.:0.00000  
##  Median :10.917   Median :46.92   Median :0.00000  
##  Mean   :10.749   Mean   :47.01   Mean   :0.09844  
##  3rd Qu.:13.750   3rd Qu.:48.75   3rd Qu.:0.00000  
##  Max.   :18.583   Max.   :52.58   Max.   :1.00000
# Last, we join this combined data set with the climate data.
sp_env <- cbind(sp_env, terra::extract(x = clim, y = sp_env[,c('decimalLongitude','decimalLatitude')], cells=T) )
##  decimalLongitude decimalLatitude      occ                ID      
##  Min.   : 2.750   Min.   :40.92   Min.   :0.00000   Min.   :   1  
##  1st Qu.: 7.583   1st Qu.:45.25   1st Qu.:0.00000   1st Qu.:1044  
##  Median :10.917   Median :46.92   Median :0.00000   Median :2088  
##  Mean   :10.749   Mean   :47.01   Mean   :0.09844   Mean   :2088  
##  3rd Qu.:13.750   3rd Qu.:48.75   3rd Qu.:0.00000   3rd Qu.:3132  
##  Max.   :18.583   Max.   :52.58   Max.   :1.00000   Max.   :4175  
##  wc2.1_10m_bio_1  wc2.1_10m_bio_2  wc2.1_10m_bio_3 wc2.1_10m_bio_4
##  Min.   :-3.347   Min.   : 5.193   Min.   :24.74   Min.   :464.7  
##  1st Qu.: 7.565   1st Qu.: 8.075   1st Qu.:31.29   1st Qu.:624.5  
##  Median : 8.936   Median : 8.561   Median :32.60   Median :657.8  
##  Mean   : 8.912   Mean   : 8.593   Mean   :32.63   Mean   :666.5  
##  3rd Qu.:10.577   3rd Qu.: 9.182   3rd Qu.:33.88   3rd Qu.:706.0  
##  Max.   :16.905   Max.   :11.446   Max.   :41.01   Max.   :819.1  
##  wc2.1_10m_bio_5  wc2.1_10m_bio_6   wc2.1_10m_bio_7 wc2.1_10m_bio_8 
##  Min.   : 7.883   Min.   :-13.749   Min.   :18.97   Min.   :-6.303  
##  1st Qu.:21.865   1st Qu.: -4.502   1st Qu.:24.88   1st Qu.: 9.468  
##  Median :23.519   Median : -2.993   Median :26.09   Median :14.234  
##  Mean   :23.403   Mean   : -2.918   Mean   :26.32   Mean   :12.448  
##  3rd Qu.:25.830   3rd Qu.: -1.127   3rd Qu.:27.67   3rd Qu.:16.565  
##  Max.   :31.365   Max.   :  8.387   Max.   :32.22   Max.   :21.641  
##  wc2.1_10m_bio_9   wc2.1_10m_bio_10 wc2.1_10m_bio_11  wc2.1_10m_bio_12
##  Min.   :-9.6958   Min.   : 3.722   Min.   :-9.7354   Min.   : 217.0  
##  1st Qu.: 0.9681   1st Qu.:15.845   1st Qu.:-0.7333   1st Qu.: 715.0  
##  Median : 3.4552   Median :17.213   Median : 0.6764   Median : 853.0  
##  Mean   : 5.5607   Mean   :17.209   Mean   : 0.9180   Mean   : 916.5  
##  3rd Qu.: 6.3344   3rd Qu.:19.298   3rd Qu.: 2.6065   3rd Qu.:1083.0  
##  Max.   :24.3333   Max.   :24.333   Max.   :11.0833   Max.   :2334.0  
##  wc2.1_10m_bio_13 wc2.1_10m_bio_14 wc2.1_10m_bio_15 wc2.1_10m_bio_16
##  Min.   : 28.0    Min.   :  4.00   Min.   : 4.41    Min.   : 77     
##  1st Qu.: 83.0    1st Qu.: 36.00   1st Qu.:16.22    1st Qu.:226     
##  Median :100.0    Median : 49.00   Median :23.07    Median :269     
##  Mean   :107.5    Mean   : 50.24   Mean   :23.89    Mean   :296     
##  3rd Qu.:126.0    3rd Qu.: 61.00   3rd Qu.:30.06    3rd Qu.:345     
##  Max.   :278.0    Max.   :119.00   Max.   :61.70    Max.   :743     
##  wc2.1_10m_bio_17 wc2.1_10m_bio_18 wc2.1_10m_bio_19      cell      
##  Min.   : 25.0    Min.   : 31.0    Min.   : 36      Min.   :41864  
##  1st Qu.:126.0    1st Qu.:195.0    1st Qu.:146      1st Qu.:49692  
##  Median :164.0    Median :229.0    Median :189      Median :53471  
##  Mean   :170.7    Mean   :252.7    Mean   :199      Mean   :53292  
##  3rd Qu.:206.0    3rd Qu.:290.0    3rd Qu.:242      3rd Qu.:56898  
##  Max.   :438.0    Max.   :573.0    Max.   :526      Max.   :65823

4.1 Spatial thinning

When we prepare our distribution data for species distribution modelling, we also need to think about spatial autocorrelation. Using adjacent cells in model building can lead to problems with spatial autocorrelation. As a rule of thumb, data points should be at least 2-3 cells apart.

One way to avoid this is spatially thinning the records, for example using the package spThin (Aiello-Lammens et al. 2015). Load the package and look up the help page ?thin.


# The spThin package requires longitude/latitude coordinates, which we already have.
# Look up the help page and try to understand the function:

# thin() expects that the data.frame contains a column with the species name
sp_env$sp <- 'Sorex_alpinus'
# Remove adjacent cells of presence/background data:
xy <- thin(sp_env, lat.col='decimalLatitude',long.col='decimalLongitude',spec.col='sp', thin.par=30,reps=1, write.files=F,locs.thinned.list.return=T)
## ********************************************** 
##  Beginning Spatial Thinning.
##  Script Started at: Mon Nov  6 09:07:37 2023
## lat.long.thin.count
## 743 
##   1 
## [1] "Maximum number of records after thinning: 743"
## [1] "Number of data.frames with max records: 1"
## [1] "No files written for this run."
# Keep the coordinates with the most presence records
xy_keep <- xy[[1]]

# Thin the dataset - here, we first extract the cell numbers for the thinned coordinates and then use these to subset our data frame.
cells_thinned <- terra::cellFromXY(clim, xy_keep)
sp_thinned <- sp_env[sp_env$cell %in% cells_thinned,]

# Plot the map and data
plot(bg, col='grey90', legend=F)
points(sp_thinned[,1:2],pch=19,col=c('black','red')[as.factor(sp_thinned$occ)], cex=0.3)

Finally, don’t forget to save your data, for example by writing the final data frame to file or by saving the R object(s).


Alternative: The thin() function can take quite long and needs a lot of memory space. Alternatively, we can use the terra package to thin the data to a checkerboard pattern, i.e. remove every second cell.

# Create checkerboard SpatRaster 
r_chess <- mask(init(region_buf,'chess'), region_buf)
values(r_chess)[values(r_chess)<1] <- NA
names(r_chess) <- 'chess'

# Thinning to the checkerboard pattern
sp_thinned2 <- merge(,cell=T),sp_env,by='cell')

# Plot the map and data
plot(bg, col='grey90', legend=F)
points(sp_thinned2[,c('decimalLongitude','decimalLatitude')],pch=19,col=c('black','red')[as.factor(sp_thinned2$occ)], cex=0.3)

We can also thin the data to an even coarser grid.

# Create checkerboard SpatRaster with 2-times coarser resolution
region_buf_agg2 <- aggregate(region_buf,2)
r_chess_agg2 <- mask(init(region_buf_agg2,'chess'), region_buf_agg2)
values(r_chess_agg2)[values(r_chess_agg2)<1] <- NA

# extract coordinates of the white fields of the checkerboard
coords_chess_agg2 <-,xy=T)[,c('x','y')]

# Get the corresponding coordinates/cells at the original resolution
cells_chess <- terra::extract(region_buf, coords_chess_agg2, cell=T)[,c(1,3)]

# Thinning to coarse-scale checkerboard pattern
sp_thinned_agg2 <-merge(cells_chess, sp_env,by='cell')

# Plot the map and data
plot(region, col='grey90', legend=F)
points(sp_thinned_agg2[,c('decimalLongitude','decimalLatitude')],pch=19,col=c('black','red')[as.factor(sp_thinned2$occ)], cex=0.3)

These different thinning methods of coarse result in very different number of presences and absences retained.

# spThin function
##   0   1 
## 697 115
# checkerboard at original spatial resolution
##    0    1 
## 1884  193
# checkerboard at coarse spatial resolution
##   0   1 
## 447  48


In practical b1, you have downloaded your own GBIF data. Carry out the pseudo-absence data selection for this species, and run a GLM analysis.

  • Decide on a pseudo-absence/background data strategy for your species (in the simplest case, just follow the strategy used in the example workflow in section 4). Prepare your dataset with presences, pseudo-absences and the corresponding climate data at those locations (following the workflow in section 4)
  • Use the data and build a GLM. Remember the different steps of checking for multicollinearity, building the full model and simplifying it.
  • Assess the predictive performance of that model.
  • Make predictions to current climate. Decide whether you want to make predictions to the entire globe or to a more restricted geographic area (e.g. the continent where your species occurs)?


Aiello-Lammens, M. E., R. A. Boria, A. Radosavljevic, B. Vilela, and R. P. Anderson. 2015. spThin: An r Package for Spatial Thinning of Species Occurrence Records for Use in Ecological Niche Models.” Ecography 38 (5): 541–45.
Barbet-Massin, M., F. Jiguet, C. H. Albert, and W. Thuiller. 2012. “Selecting Pseudo-Absences for Species Distribution Models: How, Where and How Many?” Methods in Ecology and Evolution 3 (2): 327–38.
Kramer-Schadt, S., J. Niedballa, J. D. Pilgrim, B. Schroeder, J. Lindenborn, V. Reinfelder, M. Stillfried, et al. 2013. “The Importance of Correcting for Sampling Bias in MaxEnt Species Distribution Models.” Diversity and Distributions 19 (11): 1366–79.
Phillips, S. J., M. Dudik, J. Elith, C. H. Graham, A. Lehmann, J. Leathwick, and S. Ferrier. 2009. “Sample Selection Bias and Presence-Only Distribution Models: Implications for Background and Pseudo-Absence Data.” Ecological Applications 19 (1): 181–97.