SDMaps2 is a program to implement spatial modelling of species distribution and abundance data developed at the Dutch Centre for Field Ornithology. This document aims to introduce the user into its use by working through a series of increasingly complex modelling tasks.

Your first modelling run with SDMaps2

Provide input data

For modelling with SDMaps2 you need to prepare the following items:

a table with the spatially referenced observations
a directory holding covariate grid maps
shapefile(s) with features wou want to display on the maps
if necessary, a table which relates speciescodes to full speciesnames

Observations

Open the file lancol.csv in your favourite spreadsheet programme. (“lancol” stands for Lanius collurio, the red-backed shrike, which is the species that will accompany you through the majority of examples. In modelling runs, species are often represented by numerical codes, e.g. EURING codes for birds, or abbreviations to avoid the handling of bulky species names.) This file contains the input data to SDMaps2 and consists of geo-referenced information on where the observations were made along with the observations themselves and the corresponding environmental information. In the remainder of this vignette the locations are called the ‘plots’, the data thus are the ‘plotdata’. Let us check them in detail:

The first column (“Species”) specifies the species that you want to model. Here, it is only one (lancol), but you may analyse many species in one modelling run.
The second column (“plotID”) assigns a unique code to each location where observations were made.
Column “observed” contains the observations (for example, absence or presence of the species or the number of individuals counted in the plot). Absence must be coded as 0, presence may be coded as any number larger than zero.
Columns “x” and “y” are the geographic coordinates of the plots. Any coordinate reference system may be used, for example the WGS84 coordinate reference system, the Universal Transverse Mercator coordinate system or any national reference system.
All the columns to the right of the coordinates are environmental explanatory variables (climatic variables, land-use variables etc.) characterising the plots. We call these variables covariates in the rest of this document.

There are some important rules for naming the columns: “Species”, “plotID”, “observed”, “x” and “y” must be written exactly in this way. The variable names should be without spaces or hyphens (“-“), but the use of the underscore (“_”) is allowed. This also applies to the plot ID’s and the species codes used in column “Species”.

Another advice is how to deal with missing data, when, for example, for some given plots no data are available for one or more covariates. Although this do not pose problems for some modelling techniques, others depend on complete covariate information. SDMaps contains an imputation routine for these cases, but the procedure may cause considerably longer computing times. The general advice is, thus, to avoid missing data in the covariates.

Covariate grids

You must provide a grid with the covariates that will be used as predictors, covering the entire area for which maps have to be made. We call this grid the prediction grid, or shortly, the predgrid. The maps we use in this turorial are stored in directory [your working directory]/grids10km/. Check their names: Apart from their extension, they are identical to the column headers of the covariates in file lancol.csv. Consequently, they must not include spaces or hyphens neither! As you may suspect when you look at the directory name, the maps used here are on a rather coarse grid of 10 x 10 km, which guarantees that runtime of our examples will be short. The maps can be of ASCII or GeoTIFF format. The coordinate reference system must be identical to that used for the observations.

Shapefile with features

Output maps consisting of a grid only, however, may not be very appealing or even hard to read. At least the contours of the predgrid should be made available, but you may be interested to include other informative features such as province borders, rivers etc. Its coordinate reference system has to be the same as the observations and the covariate grids.

Species codes and names

While during modelling work you may not want to deal with bulky species names, you want to see use names in map titles or output filenames to make them more readable. In this case provide a textfile called speciesnames.txt that contains two comma-separated columns “code” (the code or abbreviation you are using during modelling) and “species” (the full species name you want to use in maps and output files).

Prepare data with function data2SDMaps()

Now we can begin with our first SDMaps function, called data2SDMaps(). It reads in all necessary data and converts them into an object of class “list”. This object–we call it the sdmdata–contains all data, paths and modelling information specified by the following arguments. The object can directly be used for a modelling run or can be saved for later use and loaded into the R working space. The saved sdmdata object should not, however, be moved into another directory because of the path definitions it contains! Using it from within another directory or on another computer would cause the analysis to crash (unless you know how to edit it manually). If you want to repeat the analysis on another computer or somewhere else in the directory structure of a computer, invoke data2SDMaps() again at the new position prior to do a modelling run.

All the input data must be provided in certain directories: the plotdata, the contour shapefile and the speciesnames file should reside in the working directory, the covariate grids should be saved in a a separate directory. Here, we set the directories according to our own directory structure; you will have to define your own directories.

# Set directories
# setwd('/data/Sovon/SDMaps2/SDMaps2_examples/example1')
grids.dir <- file.path(getwd(), 'grids10km')

Let us check the input data:

# check the plotdata
plotdata <- read.csv('lancol.csv')
print(dim(plotdata))
#> [1] 5493   31
print(head(plotdata[, 1:10]))
#>   Species   plotID observed      x      y     alt bio01 bio05    bio06   bio12
#> 1  lancol  AC25131        0 156960 514565 79.0000   110   274 -41.0000 547.000
#> 2  lancol AC251312        0 157360 513765 79.0000   110   274 -41.0000 547.000
#> 3  lancol AC251313        1 157760 513765 79.0000   110   274 -41.0000 547.545
#> 4  lancol AC251314        0 158160 513765 79.0000   110   274 -41.0000 548.000
#> 5  lancol AC251319        0 158160 513365 79.0000   110   274 -41.0000 548.000
#> 6  lancol AC251321        1 156960 512965 79.8261   110   274 -40.8261 547.826
# you see the five mandatory columns Species, plotID, observed and x
# and y in columns 1 to 5, to the right you see the first 5 of 26
# explanatory variables (altitude and some bioclimatic variables)

# check the table of species code and species names
specs <- read.table('speciesnames.txt', header = TRUE, sep = ' ')
print(specs)
#>            species   code
#> 1  Lanius collurio lancol
#> 2     Lanius minor lanmin
#> 3 Lanius excubitor lanexc

# check the contour file (Romania)
library(sf)
#> Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1; sf_use_s2() is TRUE
contour <- read_sf('contour.shp', as_tibble = FALSE)
contour <- st_cast(contour, 'MULTILINESTRING')
contour <- contour[names(contour)[1]]
names(contour)[1] <- 'Contour'
plot(contour)

# check the covariate grids
print(grids.dir)
#> [1] "/data/Sovon/SDMaps2/Paketproduktion/SDMaps2/vignettes/grids10km"
list.files(file.path(getwd(), 'grids10km'))
#>  [1] "alt_categ.asc" "alt.asc"       "bio01.asc"     "bio05.asc"    
#>  [5] "bio06.asc"     "bio12.asc"     "etp.asc"       "gr01.asc"     
#>  [9] "gr02.asc"      "gr03.asc"      "gr04.asc"      "gr05.asc"     
#> [13] "gr06.asc"      "gr07.asc"      "gr08.asc"      "gr09.asc"     
#> [17] "gr10.asc"      "gr11.asc"      "gr12.asc"      "gr13.asc"     
#> [21] "gr14.asc"      "gr15.asc"      "gr16.asc"      "gr17.asc"     
#> [25] "gr18.asc"      "prec.asc"      "tmean.asc"

# show one of these grids (altitude):
library(terra)
#> terra 1.7.29
alt <- rast(file.path(getwd(), 'grids10km', 'alt.asc'))
plot(alt)

Now everything is ready for the first run. Let us make the sdmdata with function data2SDMaps2:

# Provide information on the which coordinate reference system is
# used. Here, we apply the CRS Dealul Piscului 1970/Stereo 70
# (EPSG:31700)–our example comes from Romania–which corresponds to the
# string that is assigned to object crs:
crs <- '+proj=sterea +lat_0=46 +lon_0=25 +k=0.99975 +x_0=500000 +y_0=500000 +ellps=krass +units=m +no_defs'
sdmdata <- data2SDMaps2(# we call the sdmdata simply sdmdata:
                        named = 'sdmdata',
            # the output will be save here (if the
            # directory does not yet exist it will be
            # created):
                        outdir = file.path(getwd(), 'results'),
                        plot.data = 'lancol.csv', 
                        crs = crs, 
                        grids.dir = file.path(getwd(), 'grids10km'),
            # categorical covariates have to be mentioned
            # explicitly:
                        categ.covariate = 'alt_categ')
#> 
#> 
#>  ---------------------------
#>  Preparing data for SDMaps2
#>  ---------------------------
#> 
#>  Checking input...
#>      Checking arguments...Ok
#>      Reading plot data...Ok
#>      Checking column headers...Ok
#>      Checking plotIDs for blanks, hyphens and leading zeroes...Ok
#>      Checking species codes...Ok
#>  Checking input...Ok
#> 
#>  Adding USER-mapdata to plotdata...
#>      Adding alt_categ to data... Ok
#>      Adding alt to data...already in plot data. Adding only to prediction data Ok
#>      Adding bio01 to data...already in plot data. Adding only to prediction data Ok
#>      Adding bio05 to data...already in plot data. Adding only to prediction data Ok
#>      Adding bio06 to data...already in plot data. Adding only to prediction data Ok
#>      Adding bio12 to data...already in plot data. Adding only to prediction data Ok
#>      Adding etp to data...already in plot data. Adding only to prediction data Ok
#>      Adding gr01 to data...already in plot data. Adding only to prediction data Ok
#>      Adding gr02 to data...already in plot data. Adding only to prediction data Ok
#>      Adding gr03 to data...already in plot data. Adding only to prediction data Ok
#>      Adding gr04 to data...already in plot data. Adding only to prediction data Ok
#>      Adding gr05 to data...already in plot data. Adding only to prediction data Ok
#>      Adding gr06 to data...already in plot data. Adding only to prediction data Ok
#>      Adding gr07 to data...already in plot data. Adding only to prediction data Ok
#>      Adding gr08 to data...already in plot data. Adding only to prediction data Ok
#>      Adding gr09 to data...already in plot data. Adding only to prediction data Ok
#>      Adding gr10 to data...already in plot data. Adding only to prediction data Ok
#>      Adding gr11 to data...already in plot data. Adding only to prediction data Ok
#>      Adding gr12 to data...already in plot data. Adding only to prediction data Ok
#>      Adding gr13 to data...already in plot data. Adding only to prediction data Ok
#>      Adding gr14 to data...already in plot data. Adding only to prediction data Ok
#>      Adding gr15 to data...already in plot data. Adding only to prediction data Ok
#>      Adding gr16 to data...already in plot data. Adding only to prediction data Ok
#>      Adding gr17 to data...already in plot data. Adding only to prediction data Ok
#>      Adding gr18 to data...already in plot data. Adding only to prediction data Ok
#>      Adding prec to data...already in plot data. Adding only to prediction data Ok
#>      Adding tmean to data...already in plot data. Adding only to prediction data Ok
#> 
#>   27 user defined maps added successfully!
#> 
#>  --------------------------
#>  Execution successful ! ! !
#>  --------------------------

The on-screen messages show you what happens during the execution of data2SDMaps2 and are self-explaining. A comment on the lines ‘Adding [variable] to data…’: if a covariate is already included to the plotdata (as is here the case for all but one covariate) they are only used to make the predgrid. For covariates that are not included to the plotdata (as is here the case for alt_categ) their values for the plotIDs are determined by an overlay with the corresponding grid file. Here, a short explanation of the used covariates:

Covariate code	Explanation
alt_categ	altitude in three classes
alt	altitude
bio**	bioclimatic variables
etp	evapotranspiration
gr**	CORINE land cover classes
prec	annual precipitation
tmean	annual mean temperature

By the way: alt_categ has only been used to show how categorical covariates have to be dealt with.

Next, check the output directory. Within directory “results” the directory “OBSERVATIONS” has been created, holding file lancol.DAT and another directory, “Temp1”, which holds the sdmdata file sdmdata.RDATA and file lancol.csv (which is not identical with lancol.csv that was used for data input). At first sight, this might appear as an excessively complicated structure, but it guarantees that with these sdmdata a series of different model runs for selected species or selected covariates can be performed without the necessity to prepare new specific data for each run.

Now inspect the structure of object sdmdata in the R workspace:

str(sdmdata, max.level = 2)
#> List of 11
#>  $ data          :Classes 'sf' and 'data.frame': 5493 obs. of  29 variables:
#>   ..$ plotID   : chr [1:5493] "AC25131" "AC251312" "AC251313" "AC251314" ...
#>   ..$ alt      : num [1:5493] 79 79 79 79 79 ...
#>   ..$ bio01    : num [1:5493] 110 110 110 110 110 110 110 110 110 110 ...
#>   ..$ bio05    : num [1:5493] 274 274 274 274 274 274 274 274 274 274 ...
#>   ..$ bio06    : num [1:5493] -41 -41 -41 -41 -41 ...
#>   ..$ bio12    : num [1:5493] 547 547 548 548 548 ...
#>   ..$ prec     : num [1:5493] 59 59 59 59 59 59 59 59 59 59 ...
#>   ..$ tmean    : num [1:5493] 155 155 155 155 155 155 155 155 155 155 ...
#>   ..$ etp      : num [1:5493] 9 9 9 9 9 9 9 9 9 9 ...
#>   ..$ gr01     : num [1:5493] 0 0 0 0 0 0 0 0 0 0 ...
#>   ..$ gr02     : num [1:5493] 0 0 0 0 0 0 0 0 0 0 ...
#>   ..$ gr03     : num [1:5493] 0 0 0 0 0 0 0 0 0 0 ...
#>   ..$ gr04     : num [1:5493] 1 0.5766 0.7273 0.6072 0.0195 ...
#>   ..$ gr05     : num [1:5493] 0 0 0 0 0 0 0 0 0 0 ...
#>   ..$ gr06     : num [1:5493] 0 0 0 0 0 0 0 0 0 0 ...
#>   ..$ gr07     : num [1:5493] 0 0 0 0 0 0 0 0 0 0 ...
#>   ..$ gr08     : num [1:5493] 0 0 0 0 0 0 0 0 0 0 ...
#>   ..$ gr09     : num [1:5493] 0 0 0 0 0.719 ...
#>   ..$ gr10     : num [1:5493] 0 0 0 0 0 0 0 0 0 0 ...
#>   ..$ gr11     : num [1:5493] 0 0 0 0 0 0 0 0 0 0 ...
#>   ..$ gr12     : num [1:5493] 0 0 0 0 0 0 0 0 0 0 ...
#>   ..$ gr13     : num [1:5493] 0 0 0 0 0 0 0 0 0 0 ...
#>   ..$ gr14     : num [1:5493] 0 0 0 0 0 0 0 0 0 0 ...
#>   ..$ gr15     : num [1:5493] 0 0 0 0 0 0 0 0 0 0 ...
#>   ..$ gr16     : num [1:5493] 0 0.423 0.273 0.393 0.261 ...
#>   ..$ gr17     : num [1:5493] 0 0 0 0 0 0 0 0 0 0 ...
#>   ..$ gr18     : num [1:5493] 0 0 0 0 0 0 0 0 0 0 ...
#>   ..$ geometry :sfc_POINT of length 5493; first list element:  'XY' num [1:2] 156960 514565
#>   ..$ alt_categ: Factor w/ 3 levels "0","500","1000": 1 1 1 1 1 1 1 1 1 1 ...
#>   ..- attr(*, "sf_column")= chr "geometry"
#>   ..- attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA NA NA NA NA NA NA ...
#>   .. ..- attr(*, "names")= chr [1:28] "plotID" "alt" "bio01" "bio05" ...
#>  $ pred.data     :Classes 'sf' and 'data.frame': 4800 obs. of  28 variables:
#>   ..$ alt_categ: Factor w/ 4 levels "0","500","1000",..: 4 4 4 4 4 4 4 4 4 4 ...
#>   ..$ alt      : num [1:4800] NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ...
#>   ..$ bio01    : num [1:4800] NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ...
#>   ..$ bio05    : num [1:4800] NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ...
#>   ..$ bio06    : num [1:4800] NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ...
#>   ..$ bio12    : num [1:4800] NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ...
#>   ..$ etp      : num [1:4800] NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ...
#>   ..$ gr01     : num [1:4800] NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ...
#>   ..$ gr02     : num [1:4800] NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ...
#>   ..$ gr03     : num [1:4800] NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ...
#>   ..$ gr04     : num [1:4800] NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ...
#>   ..$ gr05     : num [1:4800] NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ...
#>   ..$ gr06     : num [1:4800] NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ...
#>   ..$ gr07     : num [1:4800] NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ...
#>   ..$ gr08     : num [1:4800] NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ...
#>   ..$ gr09     : num [1:4800] NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ...
#>   ..$ gr10     : num [1:4800] NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ...
#>   ..$ gr11     : num [1:4800] NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ...
#>   ..$ gr12     : num [1:4800] NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ...
#>   ..$ gr13     : num [1:4800] NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ...
#>   ..$ gr14     : num [1:4800] NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ...
#>   ..$ gr15     : num [1:4800] NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ...
#>   ..$ gr16     : num [1:4800] NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ...
#>   ..$ gr17     : num [1:4800] NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ...
#>   ..$ gr18     : num [1:4800] NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ...
#>   ..$ prec     : num [1:4800] NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ...
#>   ..$ tmean    : num [1:4800] NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ...
#>   ..$ geometry :sfc_POINT of length 4800; first list element:  'XY' num [1:2] 105760 795765
#>   ..- attr(*, "sf_column")= chr "geometry"
#>   ..- attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA NA NA NA NA NA NA ...
#>   .. ..- attr(*, "names")= chr [1:27] "alt_categ" "alt" "bio01" "bio05" ...
#>  $ data.name     : chr "sdmdata"
#>  $ files         : chr "/data/Sovon/SDMaps2/Paketproduktion/SDMaps2/vignettes/results/OBSERVATIONS/Temp1"
#>  $ Userdir       : chr "/data/Sovon/SDMaps2/Paketproduktion/SDMaps2/vignettes/grids10km"
#>  $ n.unique.plots: int 5493
#>  $ specnames     : chr "lancol.DAT"
#>  $ not_run       : logi FALSE
#>  $ not_run_errors: logi FALSE
#>  $ call          : language data2SDMaps2(plot.data = "lancol.csv", crs = crs, outdir = file.path(getwd(),      "results"), grids.dir = file.p| __truncated__ ...
#>  $ attrs         :List of 1
#>   ..$ generated.zeroes: logi FALSE
#>  - attr(*, "SDMdata")= logi TRUE

As reported, it is a list consisting of 11 elements. The first two elements ($ data and $ pred.data) are simple features with geometry POINT, a class representing spatial data (Pebesma 2018), of the plot data and the maps used for prediction. The other elements save names, paths, the call with which data2SDMaps2 was invoked and some additional data information. With sdmdata in the workspace we can immediately proceed to a model run. However, it is wise to explore the data prior to modelling. We will do this in the next chapter.

Explore the data

Data exploration is a crucial step in every modelling attempt, because rubbish input produces rubbish output. Outliers, heterogeneity of variance, collinearity, dependence of observations etc. may cause serious problems dependent on the chosen modelling technique. In SDMmaps2, thus, function dataExploration2 checks for a series of potential violations of the underlying assumptions relevant for the employed statistical procedure, based on the protocol suggested by Zuur, Ieno, and Elphick (2010). Run function dataExploration2, below we will discuss it line by line:

dataExploration2(SDMdata = sdmdata,
                 spec.subs = c('lancol'),
                 outdir = file.path(getwd(), 'results',
                            'DATA_EXPLORATION'),
                 data.type = 'presence',
                 path.to.contour = file.path(getwd(), 'contour.shp'),
                 path.to.speciesnames = getwd(),
                 plot.location = TRUE,
                 covar.plots = TRUE,
                 pairwise.covars = TRUE,
                 variance.inflation.factors = TRUE,
                 obs.vs.covars = TRUE,
                 outlier.detection = TRUE,
                 min.corr = 0.8,
                 dec = ','
                 )
#> Make plots of covariates...
#>                              Ok
#> Make locations plot...
#>                          Ok
#> Determine correlation between covariates...Ok
#> Calculate variance inflation factors...Ok
#> Show relationhips between covariates and observed values...
#>                                                              Ok
#> Outlier detection...
#>                      Ok
#> 
#>  Data exploration stored in /data/Sovon/SDMaps2/Paketproduktion/SDMaps2/vignettes/results/DATA_EXPLORATION

With argument SDMdata (line 1) you define the sdmdata to be explored. You know already that these data must be an object generated with function data2SDMaps2.
With argument spec.subs (line 2) you specify for which species you want to make the data exploration. This is necessary because you have perhaps a large number of species represented in the sdmdata and you want to explore the data only for one or a few selected species. Here, we have only one species in the sdmdata, nevertheless you have to make the specification. The species code used in the example is lancol.
With argument outdir (line 3) you determine the path to the directory for saving the output of function dataExploration2. If the specified directory does not exist yet, it will be generated. Here, we choose to generate a subdirectory within our output directory “results” generated in data2SDMaps2 and to call it “DATA_EXPLORATION”.
Later on we will make an SDM based on presence-absence data. Visualisation of data differs if we use different data types (counts, densities, indices), so we have to specify this with argument data.type in line 4.
With path.to.contour and path.to.speciesnames (lines 5 and 6) we provide the information where to find the shapefile with the feautures such as province limits, rivers etc (here: only the contour of Romania) and the table with species codes and full species names.
From line 7 to 13 we determine which property of the data we want to explore. If you want to do all exploration steps then you can skip them and simply use argument “all = TRUE” instead. Here, we want to do all steps, too, but we will explain each argument explicitly.
Set argument plot.location (line 8) to TRUE, if you want make a map of the plots with presence and absence data. If set to TRUE do not forget to specify the path to the shape-file with the contour including its file name (argument path.to.contour in line 5). If no path is specified, SDMaps looks for a file named “contour.shp” in the working directory. If this file is not available neither, then a warning is displayed and the location plot cannot be drawn. Attention: the shp-file must be accompanied by shx- and dbf-files, e.g. contour.shp, contour.shx and contour.dbf! In this example the file “Locations_lancol_Lanius collurio.png” is produced.
With argument covar.plots (line 9) you decide whether plots of the covariates should be made. The resulting plots show the distributions of each covariate for the plot data (your plots in space with presence and absence information of lancol) and the prediction grid. You remember: the prediction grid is the complete grid–here, of Romania–for which we want to make predictions of occurrence probability). The name of the output file is “covariate_plots.pdf”.
In line 10 you decide whether pairwise covariate plots should be made and correlation coefficients calculated. When the sdmdata contain ten covariates or less, plots and correlation coefficients are represented in a png-file. For more than ten covariates–as is the case here–a csv-file will be saved instead, containing the correlation matrix (“Pairwise_covariates_lancol.csv”). For large numbers of covariates this file can be very bulky. In this case dataExploration2 also stores a table called Correlated_covariates_[species code].csv containing only covariates with an absolute correlation coefficient of equal to or larger than the value specified in argument min.corr (line 14). The default of min.corr is 0.75.
In line 11 you specifiy if variance inflation factors (VIFs) should be calculated. VIFs are a smarter way to identify variables with collinearity problems than pure correlation factors. A VIF is calculated as the square root of $1/(1–R^2)$ where $R^2$ is from a linear regression model in which one covariate is used as a response variable, and all other covariates as explanatory variables. A high $R^2$–and thus a hight VIF– in such a model means that most of the variation in the covariate is explained by all other covariates, which means that there is collinearity. There are different opinions in the literature which values should be regarded as ‘high’. Miles and Shevlin (2001) use a threshold of 4, Wooldridge and Stewart (2015) of 10; Zuur et al. (2009) suggest a value of 5 or even 3. Here, VIFs are calculated for all covariates and the covariate with the highest VIF (if >3) is removed. The procedure is iteratively repeated until no variables with a VIF>3 are left. The result is saved in file “Variance_Inflation_Factor_lancol.txt”.
Argument obs.vs.covars (line 12) lets you decide if you want to plot the observed values (here, presence and absence) against all single covariates. Currently, two files are produced. one with package ggplot and an added smoother (“Observed_vs_covars_lancol.pdf”), the other one with R standard graphics (“Observed_vs_covars_simple_lancol.pdf”).
For the detection of outliers you can use argument outlier.detection (line 13). Its output is a Cleveland dotplot.
Tables stored by dataExploration2 can conveniently be inspected with a spread sheet. In line 15 you specify if dots or commas are used as decimal symbols.

In the output directory you find the following files (in chronological order produced during the data exploration process):

In covariate_plots.pdf you find one plot for each covariate displaying its smoothed distribution in the plot data (‘Observed’) and in the prediction grid (‘Prediction’). You can use these plots to check if the plot data are a random subset of the prediction grid or if the plots show a bias with regard to certain covariates.
Locations_lancol_Lanius collurio.png is a figure of the locations of the plot data, with absences in grey and presences in red.
The full matrix of of pairwise correlation coefficients can be inspected in file Pairwise_covariates_lancol.csv. In Correlated_covars_lancol.csv you see only those covariates that show an absolute correlation coefficient of larger than 0.8 (specified by argument min.corr):

pairwise <- read.csv2(file.path(getwd(), 'results', 'DATA_EXPLORATION',
                                'Correlated_covariates_lancol.csv'))
print(pairwise)
#>       X    alt  bio01  bio05 bio06  bio12   prec  tmean    etp
#> 1   alt     NA -0.968 -0.976    NA  0.894  0.908 -0.990 -0.894
#> 2 bio01 -0.968     NA  0.978 0.871 -0.844 -0.885  0.984  0.898
#> 3 bio05 -0.976  0.978     NA    NA -0.836 -0.851  0.994  0.899
#> 4 bio06     NA  0.871     NA    NA     NA     NA     NA     NA
#> 5 bio12  0.894 -0.844 -0.836    NA     NA  0.979 -0.860     NA
#> 6  prec  0.908 -0.885 -0.851    NA  0.979     NA -0.878     NA
#> 7 tmean -0.990  0.984  0.994    NA -0.860 -0.878     NA  0.903
#> 8   etp -0.894  0.898  0.899    NA     NA     NA  0.903     NA

In the upper part of text file Variance_Inflation_Factor_lancol.txt you see the stepwise elimination of numerical covariates. Column GVIF shows the VIFs of the single covariates, below each table the covariate with the larges VIF is shown. You thus can see the entire iterative process that leads to the final set of covariates listed under “Variables meeting VIF criterium”. You can choose to reduce the number of covariates in the analysis and use only these covariates in the modelling step in order to make the analysis clearer. The modelling techniques implemented in SDMaps2, however, are not sensitive to collinearity anyway.

vif <- readLines(file.path(getwd(), 'results', 'DATA_EXPLORATION',
                           'Variance_Inflation_Factor_lancol.txt'))
print(vif[255:257])
#> [1] "Variables meeting VIF criterium:"                                                              
#> [2] " [1] \"bio06\" \"gr01\"  \"gr02\"  \"gr03\"  \"gr04\"  \"gr05\"  \"gr06\"  \"gr07\"  \"gr08\" "
#> [3] "[10] \"gr09\"  \"gr11\"  \"gr12\"  \"gr13\"  \"gr14\"  \"gr15\"  \"gr16\"  \"gr17\"  \"gr18\" "

Observed_vs_covars_simple_lancol.pdf shows all pairwise scatterplots between the observed value (presence-absence, counts, density, …) and the covariates. Observed_vs_covars_lancol.pdf contains the same information with a nicer display and, if possible, with an added smoother along with its confidence interval in each plot.
Finally, Outlier_detection_lancol.pdf allows for the identification of data points with suspect values. Covariate etp, for example, shows one data point with value 0, whereas all othe data have values of 5 and larger (see page 2 of the output in Figure Outlier_detection_lancol.pdf). You could check out if this is a measurement of data capture error and decide to eliminate it from the data.

Make an SDM with SDMaps2

SDMaps2 is the central function of package SDMaps2. It also has a large number of arguments which you can use to specify a large range of different modelling runs. Here, however, we will begin with the most basic analyses, needing only a very few arguments. Below we will explain the call line by line.

sdm.ranger.presence <- SDMaps2(
     SDMdata = sdmdata,
     model.type = 'ranger',
     ranger.bias.correction = TRUE,
     data.type = 'presence',
     out.dir = file.path(getwd(), 'results', 'ranger', 'presence_absence'),
     vars.subs = c('-etp'),
     save.text.files = 'all',
     path.to.contour = file.path(getwd(), 'contour.shp'),
     path.to.speciesnames = getwd()
     )

In line 1 we invoke the function and save its output to object sdm.ranger.presence in the workspace. You can run SDMaps2 without saving a workspace object but this will be useful for the SDMapsSummary2 command below (section.

In line 2 we tell SDMaps2 which sdmdata to use (here, object sdmdata in the R workspace).

In line 3 we define the model type. Here, we specify that the SDM must be made with the random forest algorithm ranger. Actually, this is not necessary because ranger is the default modelling technique anyway, but we state it for clarity.

Line 4: Although regarded as one of the best available machine learning methods, random forests results are biased and tend to overestimate low values and to underestimate high values (Zhang and Lu 2012; Song 2015). By default this argument is set to TRUE anyway, for clarity we explicitly set it to TRUE here. Actually the correction has been implemented to correct count models, but it also works fine for presence models (as you will see soon).

In line 5 we specify that the data should be used as presence-absence data. In this case, all values > 0 are changed to 1.

The output directory is defined in line 6: we want a directory “ranger” to be generated within the output directory already defined for function data2SDMaps2; directory “ranger” in turn contains a directory “presence_absence”.

In line 7 we skip a variable from the analysis with argument vars.subs. This argument waits for a character vector of covariate names that should be used or skipped in the modelling run. For example, vars.subs = c(“alt”, ”bio01”) includes only covariates alt and bio01; vars.subs = c(“-gr10”, “-gr11”) uses all covariates except gr10 and gr11. You cannot mix using and skipping mode of vars.subs! For example, vars.subs = c(“bio01”, ”-gr11”) makes no sense and would cause the modelling run to stop with an error message! If left unspecified, all covariates will be included. The predictor variables defined in vars.subs will be used for all species in the modelling run.

SDMaps2 saves a lot of information in the output directory and in various formats, among others as text files. This makes it easier to use the output in other programmes but also can lead to rather bulky output. That is why by default no text files are saved. Here, we determine in line 8 that all text files be saved in the output directory.

In lines 9 and 10 we provide information where to find the contour file and the speciesnames file. You encountered these arguments already in function dataExploration.

You see the screen messages that report that SDMaps2 first checks the data, makes a copy of the analysis object in the output directory (handy if your modelling run includes more than one species and crashes; the species already run until the crash do not have to be modelled again), stores the sdmdata in a filehash database, proceeds to running the model, applies the bias correction procedure, and finally produces the maps, followed (hopefully) be a success message and the information where the results have been stored. Go to the results directory and check its contents!

First, you will find a directory called “RESULTS” in directory “ranger/presence_absence”. When you run SDMaps2 again–perhaps you want to use another selection of covariates–and do not change the specifications of argument out.dir, then its results will be stored in a new directory called “RESULTS1”. If you do repeated modelling runs, directories “RESULTS2”, “RESULTS3”, and so on, will be created. Anyhow SDMaps2 informs you where the ouput has been stored and displays a corresponding message.

Now look at the contents of “RESULTS”! First notice the naming conventions used for all the files. Since we provided the textfile speciesnames.txt, the species code could be complemented by the full (here, the scientific) species name. Files that store modelling information also include the used model type (here, ranger); maps include information what type of data they display (predicted, predictedr, int.resids). Let us go through them file by file, here listed in chronological order of creation during the modelling process.

Locations_lancol_Lanius collurio.png: The first file generated is the figure of the locations of the observations which you already know from the data exploration. In case you do not run a data exploration–which we do not recommend!–the locations of the plot data are thus shown anyway.

Locations of observed presences (red) and absences (grey)

ranger_presence_lancol_Lanius collurio.RData: This file contains the model object. Its specific structure depends on the used model type, so one should consult specific documentation of the used R-commands or packages (ranger, xgboost).
lancol_obs_fitted_gam.pdf: In this file you see the fitted and observed values before bias correction (top left), after bias correction (bottom left) and the function that maps the fitted to the corrected values (bottom right).

Effect of ranger bias correction

PointPredictions_ranger_presence_lancol_Lanius collurio.txt: This file stores the plot data (plotID, observed, x-coordinates, y-coordinates, covariates) along with the model predictions and the residuals for each plot and is only saved when the argument save.text.files is set accordingly.

pointpreds <- read.table(file.path(getwd(), 'results', 'ranger',
                            'presence_absence', 'RESULTS',
                            'PointPredictions_ranger_presence_lancol_Lanius collurio.txt'),
             header = TRUE)
print(head(pointpreds, 3))
#>     plotID observed alt bio01 bio05 bio06   bio12 prec tmean etp gr01 gr02 gr03
#> 1  AC25131        0  79   110   274   -41 547.000   59   155   9    0    0    0
#> 2 AC251312        0  79   110   274   -41 547.000   59   155   9    0    0    0
#> 3 AC251313        1  79   110   274   -41 547.545   59   155   9    0    0    0
#>      gr04 gr05 gr06 gr07 gr08 gr09 gr10 gr11 gr12 gr13 gr14 gr15    gr16 gr17
#> 1 1.00000    0    0    0    0    0    0    0    0    0    0    0 0.00000    0
#> 2 0.57661    0    0    0    0    0    0    0    0    0    0    0 0.42339    0
#> 3 0.72727    0    0    0    0    0    0    0    0    0    0    0 0.27273    0
#>   gr18 alt_categ      x      y       preds       resids Species
#> 1    0         0 156960 514565 0.017173955 -0.017173955  lancol
#> 2    0         0 157360 513765 0.005329119 -0.005329119  lancol
#> 3    0         0 157760 513765 0.625261510  0.374738490  lancol

ranger_presence_lancol_raw_vs_corrected_raw_predictions.pdf: The correcting function determined with observed and fitted values (see file lancol_obs_fitted_gam.pdf) is now applied to the predictions for the entire area.

Effect of ranger bias correction

preds_ranger_presence_lancol_Lanius collurio.RData: This file contains the prediction grid which consists of a dataframe with (1) the prediction grid coordinates, (2) the values of the covariates, (3) the raw predictions (‘predictedr’) (= predictions without spatial interpolation of residuals), (4) the interpolated residuals (‘int.resids’) and (5) the predictions (‘predicted’) ( = raw predictions + interpolated residuals) as well as a column indicating the species.
int.resids_presence_lancol.asc, predictedr_presence_lancol.asc, ‘predictedr_presence_lancol_Lanius collurio.png’, predicted_presence_lancol.asc, ‘predicted_presence_lancol_Lanius collurio.png’: But what are raw predictions, predictions and interpolated residuals, that were just mentioned? Briefly, the interpolation of residuals is an approach to handle the spatial autocorrelation of the data, which is an important issue in making SDMs (Dormann et al. 2007). The raw predictions are the predictions by the model, the predictions are these raw predictions corrected by the interpolated residuals (Hengl et al. 2009). Here, we have not intended to do any correction for spatial autocorrelation (the topic will be dealt later) but SDMaps2 in any case produces these three types of maps (as ascii files and visualised as png files). However, the interpolated residuals consist only of zeroes and the predictions consequently are identical with the raw predictions.

Predicted probability of occurrence

Last not least, the file sdmdata_analysis.RDATA–we call it the analysis object–contains information on the data used for analysis and analysis-specific information such as the modelling call, paths etc. In fact it is the sdmdata object we made by data2SDMaps2 completed with information from the modelling run.

Make a model summary with SDMapsSummary2

Getting some maps is not the whole modelling story. Of course you need some information on the modelling success, parameter estimates etc. You achieve this by the SDMapsSummary2 command which again knows a number of arguments:

summary.ranger.presence <-
      SDMapsSummary2(SDMres = sdm.ranger.presence,
                     calc.varimp = TRUE,
                     make.varimp.plot = TRUE,
                     make.partial.plots = TRUE,
                     partial.plot.vars = 2,
                     make.observed.vs.fitted.plot = TRUE,
                     make.roc.plot = TRUE
                     )

With argument SDMres you specify the model for which you want to generate a summary. This can be an object resident in the workspace (we stored the output of command SDMaps2 in an object with the name sdm.ranger.presence). Alternatively, you can specify the path to the analysis object saved on the hard disk. In this case you use the argument

SDMres = 'results/ranger/presence_absence/RESULTS/sdmdata_analysis.RDATA'

Replace the string RESULTS by RESULTS2, RESULTS3 etc. if necessary. (You remember that for each analysis run using the same output directory a new RESULTS subdirectory will be generated.)

With the remaining arguments–their defaults are FALSE–you specify which information you want to be displayed.

calc.varimp and make.varimp.plot: Calculates relative importance of the covariates and plots them.
make.partial.plots and partial.plot.vars: Makes partial dependence plots. Briefly, a partial dependence function plots the change in the average predicted value as specified covariates vary over their marginal distribution [Goldstein et al. (2015)) and thus generates something like a ‘response curve’. Since it often would be too time-consuming to make partial dependence plots for all covariates, the argument partial.plot.vars allows for choosing selected covariates. You can either specify the covariates by name, e.g. by using `r partial.plot.vars = c('bio01', 'bio05', 'etp')` or by choosing the n most important covariates, e.g. by using `r partial.plot.vars = 5` which will generate partial dependence plots for the five top important explanatory variables.
make.observed.vs.fitted.plot: Here you specify if you want to plot the oberved values against the fitted ones.
make.roc.plot: With this argument you specifiy whether a receiver operating characteristic (ROC) curve should be plotted. A ROC curve illustrates the diagnostic ability of a binary classifier system by plotting sensitivity (= proportion of correctly predicted presences) against 1 – specificity (= 1 – proportion of correctly predicted absences). The area under this curve is regarded as a measure of the prediction accuracy. A measure of 1 indicates a perfect model; random predictions would yield a measure of 0.5.
There is another argument which we do not use here: calc.moran. This allows you to test the spatial autocorrelation of the model residuals by calculation Moran’s I (Bivand, Pebesma, and Gómez-Rubio 2013, 276–84). This is a considerably time-consuming process. Even for the comparably small dataset used her it costs some time, so be careful by setting this argument to TRUE!

In the screen messages you see which calculations are being made and ist last line informs you where all the products have been saved. In chronological order of generation these are:

ranger_model_summary_lancol_Lanius collurio.txt: This is the information displayed by the summary-command for an object of class ‘ranger’.
ranger_ROC-plot_lancol_Lanius collurio.pdf: In this file you can inspect the ROC curve. Models with a curve close to the top-left corner indicate a good performance. A random model is expected to give points lying along the diagonal. The area under the curve (AUC) is 0.96. The shape of the curve and the AUC indicate a good model. The model predictions for the presence of Lanius collurio range are expressed as a probability and thus are somewhere in the range between 0 and 1. What would be the optimal threshold to make a purely binary decision, present or absent? Among the various methods available to optimize the threshold there are three included in the plot: (1) the observed prevalence (proportion of plots with the species present) (‘ObsPrev’). (2) the maximixed percentage of correctly predicted presences and absences (‘MaxPCC’); (3) the maximized value of (sensitivity + specificity)/2 (‘MaxSens+Spec’), also known as the True Skill Statistic (TSS) (Allouche, Tsoar, and Kadmon 2006). It is obvious that depending on the optimization criterium the threshold values–and the resulting binary maps–can be considerably different. Among the three methods the TSS is the most recommendable, although it might perform less good when prevalence is low (Leroy et al. 2018).

ranger_Variable_Importance_lancol_Lanius collurio.csv and pdf: Random forests have an appealing way of displaying variable importance. The algorithm estimates the importance by looking at how much prediction error increases when (OOB) data for a given variable is permuted while all others are left unchanged (Liaw and Wiener 2002). The necessary calculations are carried out tree by tree during forest construction. SDMaps2 does not display the original measure of variable importance but translates it into a relative measure. Here, covariate alt (= altitude) seems to be the most important factor by far, with gr07 and gr10 being second and third. The pdf is restricted to maximal 30 covariates to keep it clearly represented, whereas the csv-file always includes all covariates used in an analysis.

ranger_Observed_vs_Fitted_lancol_Lanius collurio.pdf: This file contains three plots. First, a barplot showing the numbers of presences and absences; second, a histogram of predicted probabilities with each bar subdivided by observed values and including several types of optimal thresholds along each plot (the same you have encountered in ‘glm_ROC-plot_lancol_Lanius collurio.pdf’); and third, a goodness-of-fit plot for presence-absence data the plots are grouped into bins based on their predicted values, and the bin prevalence is calculated for each bin. For each bin also the confidence interval is plotted and the total number of plots is labeled. For more information how to interpret this plot consult the description for function calibration.plot in package PresenceAbsence (Freeman and Moisen 2008).

ranger_Partial_Dependence_Plots_lancol_Lanius collurio.pdf: Here you see the partial plots of the two most important caovariates, alt (altitude) and bio06. In general, a partial plot can show whether the relationship between a dependent variable and a covariate is linear, monotonic or more complex. If the covariates are uncorrelated, then interpretation is clear: The partial dependence plot shows how the average prediction in your dataset changes when the covariate is changed. Interpretation is, however, more complicated when the covariates are correlated. In the left figure in the pdf you see that the occurrence probability decreases with increasing altitude. The blue rugs on the x-axis show the distribution of the covariate in steps op 10 % quantiles: There are considerably more data available with low altitude and just a few high in the mountains. Neglecting the covariate distribution can be misleading by overinterpreting regions with almost no data.

Partial dependence plots

ranger_statistics.csv: In this file you find a table with the area under curve (AUC) statistics and the explained deviance. Here, we modelled only one species. When you apply an SDMaps-run to several species, the table contains the statistics for all species in the analysis, thus it wears no species code in the filename. You can open this file in a text editor, but perhaps it is more convenient to open it in a spreadsheet programme when you have several species in the analysis (that is why you can specifiy the decimal symbol in the call). The AUC you already encountered earlier in this section (AUC.sd is its standard deviation). The optimal threshold is the threshold as calculated by the TSS which you also encountered earlier when commenting the ROC-plot.

stats <- read.csv(file.path(getwd(), 'results', 'ranger',
                            'presence_absence', 'RESULTS',
                            'ranger_statistics.csv'))
print(stats)
#>                spec.code   AUC AUC.sd opt.threshold explVar
#> 1 lancol_Lanius collurio 0.998      0           0.1    30.9

SDMaps2

Christian Kampichler

2023-09-19

Introduction