R package for computing niche models using a variety of truly presence-only methods (i.e., methods that use only environmental info from the presence sites, with no background or pseudo-absence data):
- Bioclim (rectangular environmental envelope)
- Convex hull (in environmental space)
- Domain
- Mahalanobis distance
- Kernel density estimate
- Multivariate normal distribution
Barbosa A.M. & Tarroso P. (2026) proniche – an R package for presence-only niche modelling in R. Biodiversity Informatics, in press
# devtools::install_github("https://github.com/ptarroso/proniche")
# or:
# remotes::install_github("https://github.com/ptarroso/proniche")
# install.packages("terra")
# install.packages("geodata")
library(proniche)
library(terra)
library(geodata)# you can use the example data bundled with the 'proniche' package:
vars <- terra::rast(system.file("extdata/vars.tif", package = "proniche"))
occs <- read.csv(system.file("extdata/occs.csv", package = "proniche"))
# or you can download and select your own data of interest, e.g.:
vars <- geodata::worldclim_global(var = "bio", res = 5, path = "outputs")
names(vars) <- sub("wc2.1_5m_", "", names(vars))
vars <- vars[[c("bio_1", "bio_6", "bio_12", "bio_14")]]
# crop to region of interest:
vars <- terra::crop(vars, terra::ext(-10, 4, 36, 44)) # xmin, xmax, ymin, ymax
terra::plot(vars)
# import species occurrence records in this region:
occs <- geodata::sp_occurrence(genus = "Tuta",
species = "absoluta",
ext = vars,
args = "occurrenceStatus=Present")
# note that data cleaning should normally be performed at this point!
# e.g. with fuzzySim::cleanCoords() or CoordinateCleaner::clean_coordinates()
# we omit this step here, as it may be complex and it exceeds our scope
occs <- subset(occs, select = c("lon", "lat"))
# map the occurrence records:
terra::plot(vars[[1]] * 0, col = "tan", background = "lightblue",
legend = FALSE, main = "Presences")
points(occs, pch = 20, cex = 0.2)
vals <- terra::extract(vars, occs, ID = FALSE)
# remove duplicates:
nrow(vals)
#> [1] 383
vals <- unique(vals)
nrow(vals)
#> [1] 80
head(vals)
#> bio_1 bio_6 bio_12 bio_14
#> 1 14.71175 1.984000 386 4
#> 3 16.73919 5.720833 570 14
#> 7 17.62214 5.912360 268 1
#> 12 17.96979 6.000000 210 1
#> 16 16.66498 8.770213 495 1
#> 19 17.76883 6.351948 294 1par(mfrow = c(2, 2), mar = c(2, 2, 2, 1))
proniche::freqPlot(vals, vars)
bc_fit <- proniche::promodel(vals, method = "bioclim")
ch_fit <- proniche::promodel(vals, method = "convexhull")
dm_fit <- proniche::promodel(vals, method = "domain")
mm_fit <- proniche::promodel(vals, method = "mahalanobis")
km_fit <- proniche::promodel(vals, method = "kernel")
mv_fit <- proniche::promodel(vals, method = "mvnormal")The plots show the bidimensional environmental space defined by two
variables at a time. By default, the first two variables in the model
are used, but the cols argument can define a different two-variable
combination:
par(mfrow = c(2, 3))
plot(bc_fit, main = "Bioclim")
plot(ch_fit, main = "Convex Hull")
plot(dm_fit, main = "Domain")
plot(mm_fit, main = "Mahalanobis")
plot(km_fit, main = "Kernel")
plot(mv_fit, main = "Multivariate Normal")
plot(bc_fit, cols = c(2, 4), main = "Bioclim")
plot(ch_fit, cols = c(2, 4), main = "Convex Hull")
plot(dm_fit, cols = c(2, 4), main = "Domain")
plot(mm_fit, cols = c(2, 4), main = "Mahalanobis")
plot(km_fit, cols = c(2, 4), main = "Kernel")
plot(mv_fit, cols = c(2, 4), main = "Multivariate Normal")
bc_pred <- predict(bc_fit, vars)
ch_pred <- predict(ch_fit, vars)
dm_pred <- predict(dm_fit, vars)
mm_pred <- predict(mm_fit, vars)
km_pred <- predict(km_fit, vars)
mv_pred <- predict(mv_fit, vars)
par(mfrow = c(2, 3))
terra::plot(bc_pred, type = "continuous", main = "Bioclim")
terra::plot(ch_pred, type = "continuous", main = "Convex Hull")
terra::plot(dm_pred, type = "continuous", main = "Domain")
terra::plot(mm_pred, type = "continuous", main = "Mahalanobis")
terra::plot(km_pred, type = "continuous", main = "Kernel")
terra::plot(mv_pred, type = "continuous", main = "Multivariate Normal")
From Formoso-Freire et al. (2021):
“This function rescales the predictions into a ranked suitability value for each grid cell, thus informing about the suitability of climatic conditions in relative terms and hence which locations are more suitable than others. Current implementations of BIOCLIM compute a percentile distribution of the values of each environmental variable at known presence locations. Then, the closer to the 50th percentile (the median), the more suitable each location is according to that variable. However, the more variables are included in the model, the less suitable any location is considered, because it is less likely to be close to the median for all variables, and more likely to be outside the range observed at known presences for at least one variable. Our proposed rescaling of Bioclim predictions removes the dependence of the potential distribution on the number of variables included, and it provides more realistic predictions”.
bc_rcl <- proniche::quantReclass(bc_pred[[1]])
ch_rcl <- proniche::quantReclass(ch_pred[[1]])
dm_rcl <- proniche::quantReclass(dm_pred[[1]])
mm_rcl <- proniche::quantReclass(mm_pred[[1]])
km_rcl <- proniche::quantReclass(km_pred[[1]])
mv_rcl <- proniche::quantReclass(mv_pred[[1]])
par(mfrow = c(2, 3))
terra::plot(bc_rcl, range = c(0, 1), type = "continuous", main = "Bioclim")
terra::plot(ch_rcl, range = c(0, 1), type = "continuous", main = "Convex Hull")
terra::plot(dm_rcl, range = c(0, 1), type = "continuous", main = "Domain")
terra::plot(mm_rcl, range = c(0, 1), type = "continuous", main = "Mahalanobis")
terra::plot(km_rcl, range = c(0, 1), type = "continuous", main = "Kernel")
terra::plot(mv_rcl, range = c(0, 1), type = "continuous", main = "Multivariate Normal")
It is also possible to define a hard boundary for the inferred environmental niche, predicting only the pixels inside the defined envelope:
bc_fit_trunc <- truncate(bc_fit, p_obs = 1)
ch_fit_trunc <- truncate(ch_fit, p_obs = 1)
dm_fit_trunc <- truncate(dm_fit, p_obs = 1)
mm_fit_trunc <- truncate(mm_fit, p_obs = 1)
km_fit_trunc <- truncate(km_fit, p_obs = 1)
mv_fit_trunc <- truncate(mv_fit, p_obs = 1)
bc_pred_trunc <- predict(bc_fit_trunc, vars, type = "truncated")
ch_pred_trunc <- predict(ch_fit_trunc, vars, type = "truncated")
dm_pred_trunc <- predict(dm_fit_trunc, vars, type = "truncated")
mm_pred_trunc <- predict(mm_fit_trunc, vars, type = "truncated")
km_pred_trunc <- predict(km_fit_trunc, vars, type = "truncated")
mv_pred_trunc <- predict(mv_fit_trunc, vars, type = "truncated")
plot(c(bc_pred_trunc, ch_pred_trunc, dm_pred_trunc,
mm_pred_trunc, km_pred_trunc, mv_pred_trunc),
mar = c(1.5, 1.5, 2, 2))
As shown above, truncation works well for most methods based on the proportion of known presences correctly classified. The exception is Domain, because it classifies similarity to any of the presences. Thus, every presence has a similarity value of 1, because it is equal at least to itself. It is not possible in this case to set a threshold that correctly classifies a given proportion of presences, because they will all always be correctly classified. A workaround is to define the threshold using test data, for example from a buffer around the presence points. This way we have a range of climatic similarity values in the spatial vicinity of the presences and we can choose a threshold:
buf <- buffer(terra::vect(occs), 1000)
test_data <- terra::extract(vars, buf, ID = FALSE)
dm_fit_trunc_test <- truncate(dm_fit, test_data = test_data)
dm_pred_trunc_test <- predict(dm_fit_trunc_test, vars, type = "truncated")
plot(dm_pred_trunc_test, main = "Domain truncated from test data")
There are several ways to ensemble model predictions. The example below simply uses the local mean prediction across models, showing also the variance for an assessment of uncertainty:
preds <- c(bc_rcl, ch_rcl, dm_rcl, mm_rcl, km_rcl, mv_rcl)
ens_mean <- terra::app(preds, "mean")
ens_var <- terra::app(preds, "var")
par(mfrow = c(1, 2))
terra::plot(ens_mean, range = c(0, 1), main = "Ensemble mean")
terra::plot(ens_var, range = c(0, 1), main = "Ensemble variance")
The following example code lines compute several different evaluation metrics for the ensemble mean model (this section requires loading a few more packages):
library(modEvA)
library(pROC)
library(kuenm)
library(raster)
par(mfrow = c(2, 2))
# evaluation functions in the modEvA package accept spatial data directly:
modEvA::threshMeasures(obs = occs, pred = ens_mean, thresh = "maxTSS",
measures = c("CCR", "Sensitivity", "Specificity",
"kappa", "TSS"), rm.dup = TRUE, standardize = FALSE,
main = "Threshold-based metrics")
modEvA::Boyce(obs = occs, pred = ens_mean, rm.dup.points = TRUE,
main = "Continuous Boyce Index")
modEvA::AUC(obs = occs, pred = ens_mean, rm.dup = TRUE, main = "ROC curve")
# partial AUC in the pROC package requires instead a presence/absence table:
dat <- modEvA::inputMunch(obs = occs, pred = ens_mean, rm.dup = TRUE)
str(dat)
proc <- pROC::roc(response = dat$obs, predictor = dat$pred,
partial.auc = c(0, 0.2), partial.auc.correct = TRUE)
pROC::plot.roc(proc, print.auc = TRUE, auc.polygon = TRUE,
max.auc.polygon = TRUE)
# we can also compute the partial ROC ratio under a maximum potential error:
procratio <- kuenm::kuenm_proc(occ.test = occs,
model = raster::raster(ens_mean),
threshold = 5)
procratio$pROC_summary
#> Mean_AUC_ratio_at_5% pval_pROC
#> 1.418404 0.000000