Graphical Abstract

plot(stacked) # This is the end result of the project

Introduction

This is a species of ecological and ecenomic importance in both Eastern and Southern Africa. Ecologically, it (among other grazers) maintain the savannah biome in tropical Africa; feeding on grasses and being fed on by the carnivores. Economically, it is a great tourist attraction especially during their annual migration events.

Downloading the occurrence records from GBIF

library(dismo)
# occ <- gbif('Connochaetes', 'taurinus', ext = c(-17, 52, -36, 11))
head(occ)

Cleaning the occurrence data

Data from GBIF are often not ready for use in SDM. To clean it, there are many variables one can use. I will concentrate on three columns (lon, lat, and basisOfRecord). I will only retain those that are of human observation and preserved specimen under basisOfRecord. Further, I will only consider unique records (no duplicates), remove records with NAs and add species column with 1s as values.

library(tidyverse)
occ_cln <- occ |> select(lon, lat, basisOfRecord) |> 
  filter(basisOfRecord %in% c('HUMAN_OBSERVATION', 'PRESERVED_SPECIMEN')) |> 
  select(lon, lat) |> 
  unique() |>
  drop_na()

Adding species column to the data

occ_spp <- occ_cln |> mutate(species = 1)

Setting lon and lat as coordinates

This is helpful so that we can later use them to extract predictor values from bioclimatic variables.

coordinates(occ_spp) <- ~lon+lat

Downloading bioclim variables

I will download bioclim variables from worldclim database and crop them to extent of the study area, then plot one of the 19 layers with the occurrence records on.

clim <- getData('worldclim', var = 'bio', res = 10)
ext <- c(-17, 52, -36, 11)
clim_crop <- crop(clim, ext)
plot(clim_crop[[4]])
points(occ_spp)

Extracting predictor variables

We do this so that we can check for collinearity in the predictor variables and remove thos which are quite collinear.

extract <- raster::extract(clim_crop, occ_spp)

Testing for collienarity

library(usdm)
v <- vifstep(extract)

Excluding collinear variables from the clim_crop predictors

clim_used <- exclude(clim_crop, v)
clim_used

Building the sdmData

library(sdm)
d <- sdmData(species~., train = occ_spp, predictors = clim_used, bg = list(n = 1000, method = 'gRandom'))
d

Building the sdm model


model <- sdm(species~., data = d, methods = c('rf', 'svm', 'brt'), replication = c('sub', 'boot'), test.p = 30, n = 3)
Loading required package: parallel

Predicting suitable habitat using the trained model

pred <- predict(model, clim_crop)
plot(pred[[2]])

Ensembling the model output

ens <- ensemble(model, clim_crop, setting = list(method = 'weighted', stat = 'tss', opt = 2))
plot(ens)

Creating uncertainty plot around the predictions

There are always uncertainties involved in making this kind of predictions. These are expected across methods and replications. It is therefore important to report the uncertainties so that interpretation of the predictions can be more meaningful. I use standard deviation across model outputs.

pred_sd <- calc(pred, sd)
plot(pred_sd)

Stacking the uncertainty with the ensemble to plot them together

names(ens) <- 'Habitat Suitability Scores'
names(pred_sd) <- 'Std deviation around suitability'
stacked <- stack(ens, pred_sd)
plot(stacked)

Conclusion

Even though there are higher habitat suitability scores in Southern Africa than Eastern Africa, the uncertainty around the predictions are also higher in Southern Africa than Eastern Africa.

LS0tDQp0aXRsZTogPGNlbnRlcj4iU0RNIG9mIF9Db25ub2NoYWV0ZXMgdGF1cmludXNfIGluIEVhc3Rlcm4gYW5kIFNvdXRoZXJuIEFmcmljYSI8L2NlbnRlcj4NCmF1dGhvcjogIld5Y2xpZmUgQWd1bWJhIE9sdW9jaCB3eWNsaWZlb2x1b2NoQGdtYWlsLmNvbSINCm91dHB1dDogaHRtbF9ub3RlYm9vaw0KLS0tDQoNCiMgR3JhcGhpY2FsIEFic3RyYWN0DQoNCg0KYGBge3J9DQpwbG90KHN0YWNrZWQpICMgVGhpcyBpcyB0aGUgZW5kIHJlc3VsdCBvZiB0aGUgcHJvamVjdA0KYGBgDQoNCiMgSW50cm9kdWN0aW9uDQoNClRoaXMgaXMgYSBzcGVjaWVzIG9mIGVjb2xvZ2ljYWwgYW5kIGVjZW5vbWljIGltcG9ydGFuY2UgaW4gYm90aCBFYXN0ZXJuIGFuZCBTb3V0aGVybiBBZnJpY2EuIEVjb2xvZ2ljYWxseSwgaXQgKGFtb25nIG90aGVyIGdyYXplcnMpIG1haW50YWluIHRoZSBzYXZhbm5haCBiaW9tZSBpbiB0cm9waWNhbCBBZnJpY2E7IGZlZWRpbmcgb24gZ3Jhc3NlcyBhbmQgYmVpbmcgZmVkIG9uIGJ5IHRoZSBjYXJuaXZvcmVzLiBFY29ub21pY2FsbHksIGl0IGlzIGEgZ3JlYXQgdG91cmlzdCBhdHRyYWN0aW9uIGVzcGVjaWFsbHkgZHVyaW5nIHRoZWlyIGFubnVhbCBtaWdyYXRpb24gZXZlbnRzLg0KDQojIERvd25sb2FkaW5nIHRoZSBvY2N1cnJlbmNlIHJlY29yZHMgZnJvbSBbR0JJRl0oaHR0cHM6Ly93d3cuZ2JpZi5vcmcvc3BlY2llcy8yNDQxMTA1KQ0KDQpgYGB7ciBkb3dubG9hZH0NCmxpYnJhcnkoZGlzbW8pDQojIG9jYyA8LSBnYmlmKCdDb25ub2NoYWV0ZXMnLCAndGF1cmludXMnLCBleHQgPSBjKC0xNywgNTIsIC0zNiwgMTEpKQ0KaGVhZChvY2MpDQpgYGANCg0KIyBDbGVhbmluZyB0aGUgb2NjdXJyZW5jZSBkYXRhDQoNCkRhdGEgZnJvbSBHQklGIGFyZSBvZnRlbiBub3QgcmVhZHkgZm9yIHVzZSBpbiBTRE0uIFRvIGNsZWFuIGl0LCB0aGVyZSBhcmUgbWFueSB2YXJpYWJsZXMgb25lIGNhbiB1c2UuIEkgd2lsbCBjb25jZW50cmF0ZSBvbiB0aHJlZSBjb2x1bW5zIChsb24sIGxhdCwgYW5kIGJhc2lzT2ZSZWNvcmQpLiBJIHdpbGwgb25seSByZXRhaW4gdGhvc2UgdGhhdCBhcmUgb2YgaHVtYW4gb2JzZXJ2YXRpb24gYW5kIHByZXNlcnZlZCBzcGVjaW1lbiB1bmRlciBiYXNpc09mUmVjb3JkLiBGdXJ0aGVyLCBJIHdpbGwgb25seSBjb25zaWRlciB1bmlxdWUgcmVjb3JkcyAobm8gZHVwbGljYXRlcyksIHJlbW92ZSByZWNvcmRzIHdpdGggTkFzIGFuZCBhZGQgc3BlY2llcyBjb2x1bW4gd2l0aCAxcyBhcyB2YWx1ZXMuDQoNCmBgYHtyIGNsZWFuaW5nfQ0KbGlicmFyeSh0aWR5dmVyc2UpDQpvY2NfY2xuIDwtIG9jYyB8PiBzZWxlY3QobG9uLCBsYXQsIGJhc2lzT2ZSZWNvcmQpIHw+IA0KICBmaWx0ZXIoYmFzaXNPZlJlY29yZCAlaW4lIGMoJ0hVTUFOX09CU0VSVkFUSU9OJywgJ1BSRVNFUlZFRF9TUEVDSU1FTicpKSB8PiANCiAgc2VsZWN0KGxvbiwgbGF0KSB8PiANCiAgdW5pcXVlKCkgfD4NCiAgZHJvcF9uYSgpDQpgYGANCg0KIyBBZGRpbmcgc3BlY2llcyBjb2x1bW4gdG8gdGhlIGRhdGENCg0KYGBge3Igc3BlY2llc30NCm9jY19zcHAgPC0gb2NjX2NsbiB8PiBtdXRhdGUoc3BlY2llcyA9IDEpDQpgYGANCg0KIyBTZXR0aW5nIGxvbiBhbmQgbGF0IGFzIGNvb3JkaW5hdGVzDQoNClRoaXMgaXMgaGVscGZ1bCBzbyB0aGF0IHdlIGNhbiBsYXRlciB1c2UgdGhlbSB0byBleHRyYWN0IHByZWRpY3RvciB2YWx1ZXMgZnJvbSBiaW9jbGltYXRpYyB2YXJpYWJsZXMuIA0KDQpgYGB7ciBjb29yZGluYXRlc30NCmNvb3JkaW5hdGVzKG9jY19zcHApIDwtIH5sb24rbGF0DQpgYGANCg0KIyBEb3dubG9hZGluZyBiaW9jbGltIHZhcmlhYmxlcw0KDQpJIHdpbGwgZG93bmxvYWQgYmlvY2xpbSB2YXJpYWJsZXMgZnJvbSB3b3JsZGNsaW0gZGF0YWJhc2UgYW5kIGNyb3AgdGhlbSB0byBleHRlbnQgb2YgdGhlIHN0dWR5IGFyZWEsIHRoZW4gcGxvdCBvbmUgb2YgdGhlIDE5IGxheWVycyB3aXRoIHRoZSBvY2N1cnJlbmNlIHJlY29yZHMgb24uIA0KDQpgYGB7ciBjbGltfQ0KY2xpbSA8LSBnZXREYXRhKCd3b3JsZGNsaW0nLCB2YXIgPSAnYmlvJywgcmVzID0gMTApDQpleHQgPC0gYygtMTcsIDUyLCAtMzYsIDExKQ0KY2xpbV9jcm9wIDwtIGNyb3AoY2xpbSwgZXh0KQ0KcGxvdChjbGltX2Nyb3BbWzRdXSkNCnBvaW50cyhvY2Nfc3BwKQ0KYGBgDQoNCiMgRXh0cmFjdGluZyBwcmVkaWN0b3IgdmFyaWFibGVzDQoNCldlIGRvIHRoaXMgc28gdGhhdCB3ZSBjYW4gY2hlY2sgZm9yIGNvbGxpbmVhcml0eSBpbiB0aGUgcHJlZGljdG9yIHZhcmlhYmxlcyBhbmQgcmVtb3ZlIHRob3Mgd2hpY2ggYXJlIHF1aXRlIGNvbGxpbmVhci4NCg0KYGBge3IgZXh0cmFjdH0NCmV4dHJhY3QgPC0gcmFzdGVyOjpleHRyYWN0KGNsaW1fY3JvcCwgb2NjX3NwcCkNCmBgYA0KDQojIFRlc3RpbmcgZm9yIGNvbGxpZW5hcml0eQ0KDQpgYGB7ciBjb2xsaW5lYXJ9DQpsaWJyYXJ5KHVzZG0pDQp2IDwtIHZpZnN0ZXAoZXh0cmFjdCkNCmBgYA0KDQojIEV4Y2x1ZGluZyBjb2xsaW5lYXIgdmFyaWFibGVzIGZyb20gdGhlIGNsaW1fY3JvcCBwcmVkaWN0b3JzDQoNCmBgYHtyIHJlbW92ZWR9DQpjbGltX3VzZWQgPC0gZXhjbHVkZShjbGltX2Nyb3AsIHYpDQpjbGltX3VzZWQNCmBgYA0KDQojIEJ1aWxkaW5nIHRoZSBzZG1EYXRhDQoNCmBgYHtyIHNkbURhdGF9DQpsaWJyYXJ5KHNkbSkNCmQgPC0gc2RtRGF0YShzcGVjaWVzfi4sIHRyYWluID0gb2NjX3NwcCwgcHJlZGljdG9ycyA9IGNsaW1fdXNlZCwgYmcgPSBsaXN0KG4gPSAxMDAwLCBtZXRob2QgPSAnZ1JhbmRvbScpKQ0KZA0KYGBgDQoNCiMgQnVpbGRpbmcgdGhlIHNkbSBtb2RlbA0KDQpgYGB7ciBzZG19DQoNCm1vZGVsIDwtIHNkbShzcGVjaWVzfi4sIGRhdGEgPSBkLCBtZXRob2RzID0gYygncmYnLCAnc3ZtJywgJ2JydCcpLCByZXBsaWNhdGlvbiA9IGMoJ3N1YicsICdib290JyksIHRlc3QucCA9IDMwLCBuID0gMykNCm1vZGVsDQpgYGANCg0KIyBQcmVkaWN0aW5nIHN1aXRhYmxlIGhhYml0YXQgdXNpbmcgdGhlIHRyYWluZWQgbW9kZWwNCg0KYGBge3IgcHJlZGljdGlvbn0NCnByZWQgPC0gcHJlZGljdChtb2RlbCwgY2xpbV9jcm9wKQ0KcGxvdChwcmVkW1syXV0pDQpgYGANCg0KIyBFbnNlbWJsaW5nIHRoZSBtb2RlbCBvdXRwdXQNCg0KYGBge3IgZW5zZW1ibGV9DQplbnMgPC0gZW5zZW1ibGUobW9kZWwsIGNsaW1fY3JvcCwgc2V0dGluZyA9IGxpc3QobWV0aG9kID0gJ3dlaWdodGVkJywgc3RhdCA9ICd0c3MnLCBvcHQgPSAyKSkNCnBsb3QoZW5zKQ0KYGBgDQoNCiMgQ3JlYXRpbmcgdW5jZXJ0YWludHkgcGxvdCBhcm91bmQgdGhlIHByZWRpY3Rpb25zDQpUaGVyZSBhcmUgYWx3YXlzIHVuY2VydGFpbnRpZXMgaW52b2x2ZWQgaW4gbWFraW5nIHRoaXMga2luZCBvZiBwcmVkaWN0aW9ucy4gVGhlc2UgYXJlIGV4cGVjdGVkIGFjcm9zcyBtZXRob2RzIGFuZCByZXBsaWNhdGlvbnMuIEl0IGlzIHRoZXJlZm9yZSBpbXBvcnRhbnQgdG8gcmVwb3J0IHRoZSB1bmNlcnRhaW50aWVzIHNvIHRoYXQgaW50ZXJwcmV0YXRpb24gb2YgdGhlIHByZWRpY3Rpb25zIGNhbiBiZSBtb3JlIG1lYW5pbmdmdWwuIEkgdXNlIHN0YW5kYXJkIGRldmlhdGlvbiBhY3Jvc3MgbW9kZWwgb3V0cHV0cy4NCg0KYGBge3Igc3RhbmRhcmRfZGV2aWF0aW9ufQ0KcHJlZF9zZCA8LSBjYWxjKHByZWQsIHNkKQ0KcGxvdChwcmVkX3NkKQ0KYGBgDQoNCiMgU3RhY2tpbmcgdGhlIHVuY2VydGFpbnR5IHdpdGggdGhlIGVuc2VtYmxlIHRvIHBsb3QgdGhlbSB0b2dldGhlcg0KDQpgYGB7ciBzdGFja30NCm5hbWVzKGVucykgPC0gJ0hhYml0YXQgU3VpdGFiaWxpdHkgU2NvcmVzJw0KbmFtZXMocHJlZF9zZCkgPC0gJ1N0ZCBkZXZpYXRpb24gYXJvdW5kIHN1aXRhYmlsaXR5Jw0Kc3RhY2tlZCA8LSBzdGFjayhlbnMsIHByZWRfc2QpDQpwbG90KHN0YWNrZWQpDQpgYGANCg0KIyBDb25jbHVzaW9uDQoNCkV2ZW4gdGhvdWdoIHRoZXJlIGFyZSBoaWdoZXIgaGFiaXRhdCBzdWl0YWJpbGl0eSBzY29yZXMgaW4gU291dGhlcm4gQWZyaWNhIHRoYW4gRWFzdGVybiBBZnJpY2EsIHRoZSB1bmNlcnRhaW50eSBhcm91bmQgdGhlIHByZWRpY3Rpb25zIGFyZSBhbHNvIGhpZ2hlciBpbiBTb3V0aGVybiBBZnJpY2EgdGhhbiBFYXN0ZXJuIEFmcmljYS4g