Species distribution models (SDMs) have been developed for several years to adress conservation issues, assess the direct impact of human activities on ecosystems and predict the potential distribution shifts of invasive species (see Elith et al. 2006, Elith & Leathwick 2009, Pearson 2007 for reviews). SDM relate species occurrences with environmental information and can extrapolate species distribution on their entire occupied space. Applying SDM on limited occurrence datasets can therefore bring complementary information in non-visited areas. However, users must be aware of potential bias and limitations while using such poor and historical datasets (Araujo & Guisan 2006, Robinson 2011, Proosdij et al. 2016) that require corrections (Phillips et al. 2009, Barbet-Massin et al. 2012).
SDMPlay is a pedagogic package that will allow you to compute SDMs with two popular machine learning approaches, BRT (Boosted Regression Trees) and MaxEnt (Maximum Entropy). It contains occurrences of marine species and environmental descriptors datasets as examples for a first use of SDMs. You can also upload your own dataset. Basic approaches for model calibration and execution are provided. Classic tools to evaluate model performance are supplied (Area Under the Curve, omission rate and confusion matrix) and are completed with tools to perform null models (Raes & ter Steege 2007, Proosdij et al. 2016).
The biological dataset includes original occurrences of two echinoid species (sea urchins) present on the Kerguelen Plateau. The environmental dataset compiles 15 environmental descriptors, displayed in a raster format, on the extent of the Kerguelen Plateau, for different time periods.
Remarks
This package focusses on datasets containing presence-only data. Functions must be adapted whether dealing with presence-absence or abundance data. Presence-only methods imply using background data to be selected in the study area (Pearce & Boyce 2006) in order to calibrate the model. Several background sampling methods exist (Phillips et al. 2009), and its choice depends on the presence-only sampling pattern and the scientific questions. This package focusses only on a random sampling of background data. You can refere to other packages such as biomod2 if you want to use another sampling strategy (e.g. relative envelope, distance to disk, independent strategy).
In the package, you can download occurrence data of two echinoid species, Brisaster antarcticus and Ctenocidaris nutrix, distributed on the Kerguelen Plateau.These two species present contrasting ecological niches, with different feeding preferences and reproductive behaviours (David et al. 2005). The complete dataset of Kerguelen echinoid species is available in Guillaumot et al. (2016).
library(SDMPlay)
data("ctenocidaris.nutrix")
head(ctenocidaris.nutrix)
## id scientific.name scientific.name.authorship
## 56 1 Ctenocidaris_nutrix (Thomson 1876)
## 57 2 Ctenocidaris_nutrix (Thomson 1876)
## 58 3 Ctenocidaris_nutrix (Thomson 1876)
## 59 4 Ctenocidaris_nutrix (Thomson 1876)
## 60 5 Ctenocidaris_nutrix (Thomson 1876)
## 61 6 Ctenocidaris_nutrix (Thomson 1876)
## genus family
## 56 Ctenocidaris Mortensen 1910 Ctenocidarinae Mortensen 1928
## 57 Ctenocidaris Mortensen 1910 Ctenocidarinae Mortensen 1928
## 58 Ctenocidaris Mortensen 1910 Ctenocidarinae Mortensen 1928
## 59 Ctenocidaris Mortensen 1910 Ctenocidarinae Mortensen 1928
## 60 Ctenocidaris Mortensen 1910 Ctenocidarinae Mortensen 1928
## 61 Ctenocidaris Mortensen 1910 Ctenocidarinae Mortensen 1928
## order.and.higher.taxonomic.rank decimal.Longitude decimal.Latitude depth
## 56 Cidaroida Claus 1880 67.13167 -48.98500 315
## 57 Cidaroida Claus 1880 67.33167 -49.44167 301
## 58 Cidaroida Claus 1880 67.51167 -49.00500 206
## 59 Cidaroida Claus 1880 67.54167 -48.11667 365
## 60 Cidaroida Claus 1880 67.88500 -49.46667 191
## 61 Cidaroida Claus 1880 68.05833 -49.06667 178
## year campaign reference vessel
## 56 1975 MD04 (BENTHOS) De Ridder et al. 1992 Marion Dufresne
## 57 1975 MD04 (BENTHOS) De Ridder et al. 1992 Marion Dufresne
## 58 1975 MD04 (BENTHOS) De Ridder et al. 1992 Marion Dufresne
## 59 1975 MD04 (BENTHOS) De Ridder et al. 1992 Marion Dufresne
## 60 1975 MD04 (BENTHOS) De Ridder et al. 1992 Marion Dufresne
## 61 1975 MD04 (BENTHOS) De Ridder et al. 1992 Marion Dufresne
This package also contains stacks of raster layers, corresponding to environmental descriptors in the region of the Kerguelen Plateau, for three time periods [1965-1974], [2005-2012], and for the climatic scenario A1B (IPCC, 4th report 2007) for 2200. Grid-cells are set at a 0.1° resolution and data were not interpolated (presence of N/A values in the area). Extra metadata and environmental layers are available in Guillaumot et al. (2016).
Load the raster stacks
data("predictors1965_1974")
data("predictors2005_2012")
data("predictors2200AIB")
Plot the layers and explore their properties
library(raster)
plot(subset(predictors2005_2012, c(1:4)))
As you can notice, particularly for seafloor layers, maps are incomplete and contains an important number of missing values (N/A) because data were not interpolated in space. You can interpolate your data using the functions provided in the raster package, being aware of the interpretation issues related to this interpolation.
predictors2005_2012
## class : RasterStack
## dimensions : 100, 179, 17900, 15 (nrow, ncol, ncell, nlayers)
## resolution : 0.1, 0.1 (x, y)
## extent : 63, 80.9, -56, -46 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
## names : depth, seasurfac//_2005_2012, seasurfac//_2005_2012, seafloor_//_2005_2012, seafloor_//_2005_2012, seasurfac//_2005_2012, seasurfac//_2005_2012, seafloor_//_2005_2012, seafloor_//_2005_2012, chlorophy//_2005_2012, geomorphology, sediments, slope, seafloor_//_2005_2012, roughness
## min values : -4.977000e+03, 3.263100e-01, -4.002820e+00, -2.985100e-01, -3.097540e+00, 3.362941e+01, -6.510162e-02, 3.364861e+01, -1.297913e-01, 1.358515e-01, 1.000000e+00, 2.000000e+00, 4.822916e-05, 4.008017e+00, 7.000000e+00
## max values : -1.0000000, 9.0129099, -1.2873000, 4.8823700, 0.9008052, 33.9671898, 0.1001129, 34.8030014, 0.1012001, 2.7323778, 19.0000000, 13.0000000, 0.1546878, 7.6222906, 3417.0000000
names(predictors2005_2012)
## [1] "depth"
## [2] "seasurface_temperature_mean_2005_2012"
## [3] "seasurface_temperature_amplitude_2005_2012"
## [4] "seafloor_temperature_mean_2005_2012"
## [5] "seafloor_temperature_amplitude_2005_2012"
## [6] "seasurface_salinity_mean_2005_2012"
## [7] "seasurface_salinity_amplitude_2005_2012"
## [8] "seafloor_salinity_mean_2005_2012"
## [9] "seafloor_salinity_amplitude_2005_2012"
## [10] "chlorophyla_summer_mean_2005_2012"
## [11] "geomorphology"
## [12] "sediments"
## [13] "slope"
## [14] "seafloor_oxygen_mean_2005_2012"
## [15] "roughness"
The first step after checking your data is to adapt your dataset for modelling. Model algorithms require a table containing the environmental values associated with each occurrence data.
| ID* | Longitude | Latitude | depth | … | temperature |
|---|---|---|---|---|---|
| 1 | 63.33 | -48.26 | -480 | … | 1.4 |
| 1 | 64.13 | -48.57 | -104 | … | 1.2 |
| … | … | … | … | … | … |
| 0 | 67.32 | -47.23 | -1013 | … | 2.5 |
| 0 | 67.90 | -55.45 | -98 | … | 4.3 |
*ID corresponds to presence data (ID=1), or background data (ID=0).
This table will be included afterwards within the SDM algorithm (BRT, MaxEnt) and will infere presence probabilities on the area that you will define for extrapolation.
The package will guide you to build this dataframe and save it as a SDMtab object that will be loaded in the following functions of the package.
###Extract latitude and longitude values
ctenocidaris.nutrix.occ <- ctenocidaris.nutrix[,c(7,8)]
head(ctenocidaris.nutrix.occ)
## decimal.Longitude decimal.Latitude
## 56 67.13167 -48.98500
## 57 67.33167 -49.44167
## 58 67.51167 -49.00500
## 59 67.54167 -48.11667
## 60 67.88500 -49.46667
## 61 68.05833 -49.06667
SDMtable_ctenocidaris <- SDMtab(xydata=ctenocidaris.nutrix.occ,
predictors=predictors2005_2012,
unique.data=FALSE,
same=TRUE)
unique.data indicates that the function will look for presence-only data that fall on a same grid-cell pixel. When unique.data= TRUE, these presence-only duplicates will be removed from the xydata variable. same and background.nb functions refer to the sampling of background data. background.nb, indicates the specific number of background data to sample, while same is a shortcut that induces the sampling of a number of background data similar to the number of presence-only data available. You can refer to Barbet-Massin et al. (2012) to choose the most appropriate number of background data to sample for your case study.
We can display the beginning and the end of the first columns of this new SDMtab object:
head(SDMtable_ctenocidaris[,c(1:5)])
## id longitude latitude depth seasurface_temperature_mean_2005_2012
## 1 1 67.15 -48.95 -653 4.24109
## 2 1 67.35 -49.45 -204 3.89770
## 3 1 67.55 -49.05 -168 4.06841
## 4 1 67.55 -48.15 -355 4.65109
## 5 1 67.85 -49.45 -136 3.86259
## 6 1 68.05 -49.05 -155 4.03309
tail(SDMtable_ctenocidaris[,c(1:5)])
## id longitude latitude depth seasurface_temperature_mean_2005_2012
## 245 0 63.95 -48.15 -3260 4.60371
## 246 0 76.35 -54.85 -1476 1.60959
## 247 0 80.65 -51.75 -3799 2.15615
## 248 0 64.55 -51.45 -4130 3.39549
## 249 0 72.85 -47.95 -1332 4.57299
## 250 0 63.75 -48.75 -2617 4.25055
The dataframe combines environmental values of the 125 presence-only data available (ID=1) and environmental values associated with 125 background data randomly sampled in the area (ID=0).
You can display the sampled data on a map:
# nice colors
bluepalette<-colorRampPalette(c("blue4","blue","dodgerblue", "deepskyblue","lightskyblue"))(800)
# Isolate depth layer from the environmental stack
depth <- subset(predictors2005_2012,1)
# Extract background coordinates from SDMtable
background.occ <- subset(SDMtable_ctenocidaris,SDMtable_ctenocidaris$id==0)[,c(2,3)]
# plot the result on depth layer
plot(depth, col=bluepalette, cex=0.8,legend.width=0.5, legend.shrink=0.4,
legend.args=list(text='Depth (m)', side=3, font=2, cex=0.8))
points(ctenocidaris.nutrix.occ, pch= 20, col="black")
points(background.occ, pch= 20, col="red")
legend("bottomleft", pch=20, col=c("black", "red"), legend=c("presence-only data","background data"), cex=0.6)
You can assess the quality of your dataset with the SDMdata.quality function. This function estimates the percentage of presence-only data that fall on grid-cell pixels containing non-informative values (N/A). It estimates the quality of your dataset when raster layers are not interpolated.
SDMdata.quality(SDMtable_ctenocidaris)
## NA.percent (%)
## depth 1.2
## seasurface_temperature_mean_2005_2012 6.0
## seasurface_temperature_amplitude_2005_2012 6.0
## seafloor_temperature_mean_2005_2012 31.2
## seafloor_temperature_amplitude_2005_2012 31.2
## seasurface_salinity_mean_2005_2012 6.0
## seasurface_salinity_amplitude_2005_2012 4.0
## seafloor_salinity_mean_2005_2012 31.2
## seafloor_salinity_amplitude_2005_2012 31.2
## chlorophyla_summer_mean_2005_2012 6.4
## geomorphology 3.6
## sediments 2.0
## slope 0.0
## seafloor_oxygen_mean_2005_2012 44.4
## roughness 2.0
A last calibration step that you can perform before modelling is delineating the modelled area. The delim.area function can be used to restrict in geography and depth the environmental descriptors layers. This step can play an important role to enhance modelling performances by limiting the extent of extrapolation.
par(mar=c(0,0,0,0))
# restrict to 1500m depth
predictors2005_2012_1500m <- delim.area(predictors2005_2012, longmin=62, longmax=80,latmin=-55 , latmax=-45, interval=c(0,-1500))
# plot the new layer
plot(subset(predictors2005_2012_1500m,1), col=bluepalette,legend.width=0.5, legend.shrink=0.25,
legend.args=list(text='Depth (m)', side=3, font=2, cex=0.8))
You can focus your background sampling on this restrained environment. Run again the SDMtab code with these changes. The function will omit the N/A pixels when selecting the random background data.
SDMtable_ctenocidaris_1500 <- SDMtab(xydata=ctenocidaris.nutrix.occ,
predictors=predictors2005_2012_1500m,
unique.data=FALSE,
same=TRUE)
Observe the changes
Once you have built your SDMtab dataframe, you can easily perform models using the compute.maxent and compute.brt functions.
compute.brt(x, proj.predictors, tc = 2, lr = 0.001, bf = 0.75,
n.trees = 50, step.size = n.trees)
compute.maxent(x, proj.predictors)
The fonctions require two main parameters, x which correspond to the SDMtab object previously created and proj.predictors, the RasterStack containing the environmental descriptors on which you want to project your model.
The other arguments aim at calibrating the model. You can refere to Elith et al. (2008) and Elith et al. (2011) to choose the parameters according to your dataset. BRT arguments are explained in gbm package.
Extrapolate species distribution on the Kerguelen Plateau, for [2005-2012]
Cteno_model_2005_2012 <- SDMPlay:::compute.brt(x=SDMtable_ctenocidaris_1500, proj.predictors=predictors2005_2012_1500m, tc = 2, lr = 0.001, bf = 0.75, n.trees = 500)
While the function is uploading, you can observe that the gbm function, called by SDMPlay, calculates the regression trees until reaching the best estimation. This can help you refine your model calibration. Run your model while changing the calibration until reaching the best modelling performances.
Afterwards, different outputs can be produced.
# display nice colors
palettecolor <- colorRampPalette(c("deepskyblue", "darkseagreen","lightgreen","green","yellow","gold","orange", "red","firebrick"))(100)
# plot the results
plot(Cteno_model_2005_2012$raster.prediction,col=palettecolor, main="Projection for [2005-2012]",
cex.axis= 0.7,
legend.width=0.5, legend.shrink=0.25,
legend.args=list(text='Distribution probability', side=3, font=2, cex=0.8))
The output of your model cannot extrapolate on the grid-cell pixels from which it does not know environmental values. Choose the option of interpolating your RasterStack layers before modelling or when projecting if you want to obtain smoother prediction maps.
The map gives you the species distribution probabilities contained between 0 and 1.
contributions <- Cteno_model_2005_2012$response$contributions
b <- barplot(contributions[,2], ylab="contribution (%)")
text(b-0.1, par("usr")[3] - 0.025, srt = 45, adj = 1, labels=contributions[,1],cex=0.5,xpd=T)
Response plots are useful indicators of environmental preferential values for the species. y axis contains distribution probabilities predicted by the model and associates these values with environmental data.
library(dismo)
gbm.plot(Cteno_model_2005_2012$response,n.plots=12,cex.axis=0.6,cex.lab=0.7, smooth=TRUE)
Display the interaction between your environmental variables and plot them in 3D
interactions <- gbm.interactions(Cteno_model_2005_2012$response)
head(interactions$rank.list[,c(5,2,4)])
## int.size var1.names
## 1 7.08 seasurface_temperature_amplitude_2005_2012
## 2 3.02 seasurface_salinity_mean_2005_2012
## 3 2.86 seasurface_temperature_mean_2005_2012
## 4 1.80 chlorophyla_summer_mean_2005_2012
## 5 1.76 chlorophyla_summer_mean_2005_2012
## 6 0.88 seasurface_salinity_amplitude_2005_2012
## var2.names
## 1 seasurface_temperature_mean_2005_2012
## 2 seasurface_temperature_amplitude_2005_2012
## 3 depth
## 4 seasurface_temperature_mean_2005_2012
## 5 depth
## 6 seasurface_temperature_amplitude_2005_2012
gbm.perspec(Cteno_model_2005_2012$response,interactions$rank.list[1,1],interactions$rank.list[1,3], cex.lab=0.6, cex.axis=0.6,par(mar=c(0,0,0,0)))
If you want to project your model on another time period and infere your species distribution for other environmental conditions, you just need to change the proj.predictors in compute.brt. The fonction will do the relationship between the environmental descriptors used for modelling and projecting. You must ensure that the extent, order and names of your raster layers are similar.
The procedure for MaxEnt algorithm is similar to BRT. compute.maxent uses the functionalities of the dismo maxent function. This function calls MaxEnt species distribution software, which is a java program that could be downloaded here. In order to run compute.maxent, put the maxent.jar file downloaded at this address in the java folder of the dismo package (path obtained with the system.file('java', package='dismo') command).
For issues with Java installation, consult the dismo and rJava packages.
MaxEnt model outputs are similar to BRT, you can compute maps, response plots, environmental descriptors contributions. Refere to the example section of the function for details.
SDMPlay provides extra fonctions to go further in your modelling work. You can perform null models with null.model, evaluate modelling performance and define probability threshold with SDMeval. Further reading and examples are provided within the functions, don't hesitate to explore them.
Araujo, M. B., & Guisan, A. (2006). Five (or so) challenges for species distribution modelling. Journal of biogeography, 33(10), 1677-1688.
Barbet‐Massin, M., Jiguet, F., Albert, C. H., & Thuiller, W. (2012). Selecting pseudo‐absences for species distribution models: how, where and how many? Methods in Ecology and Evolution, 3(2), 327-338.
David, B., Choné, T., Mooi, R. & de Ridder C. (2005). Antarctic echinoidea (Vol. 10). ARG Gantner.
Elith, J.,P Anderson, R., Dudík, M., Ferrier, S., Guisan, A., J Hijmans, R., Huettmann, F., … & A Loiselle, B. (2006). Novel methods improve prediction of species’ distributions from occurrence data. Ecography, 29(2), 129-151.
Elith, J., Leathwick, J. R., & Hastie, T. (2008). A working guide to boosted regression trees. Journal of Animal Ecology, 77(4), 802-813.
Elith, J., & Leathwick, J. R. (2009). Species distribution models: ecological explanation and prediction across space and time. Annual review of ecology, evolution, and systematics, 40, 677-697.
Elith, J., Phillips, S. J., Hastie, T., Dudík, M., Chee, Y. E., & Yates, C. J. (2011). A statistical explanation of MaxEnt for ecologists. Diversity and distributions, 17(1), 43-57.
Guillaumot, C., Martin, A., Fabri-Ruiz, S., Eléaume, M., & Saucède, T. (2016). Echinoids of the Kerguelen Plateau–occurrence data and environmental setting for past, present, and future species distribution modelling. ZooKeys, (630), 1.
Pearce, J. L., & Boyce, M. S. (2006). Modelling distribution and abundance with presence‐only data. Journal of applied ecology, 43(3), 405-412.
Pearson, R. G. (2007). Species’ distribution modeling for conservation educators and practitioners. Synthesis. American Museum of Natural History, 50.
Phillips, S. J., Dudík, M., Elith, J., Graham, C. H., Lehmann, A., Leathwick, J., & Ferrier, S. (2009). Sample selection bias and presence‐only distribution models: implications for background and pseudo‐absence data. Ecological applications, 19(1), 181-197.
Proosdij, A. S., Sosef, M. S., Wieringa, J. J., & Raes, N. (2016). Minimum required number of specimen records to develop accurate species distribution models. Ecography, 39(6), 542-552.
Raes, N., & ter Steege, H. (2007). A null‐model for significance testing of presence‐only species distribution models. Ecography, 30(5), 727-736.
Robinson, L. M., Elith, J., Hobday, A. J., Pearson, R. G., Kendall, B. E., Possingham, H. P., & Richardson, A. J. (2011). Pushing the limits in marine species distribution modelling: lessons from the land present challenges and opportunities. Global Ecology and Biogeography, 20(6), 789-802.