R/threshold_3D.R
threshold_3D.RdCreates 3D presence/absence layers by setting all values in a suitability layer above a certain threshold to 1 and all values below that threshold to 0. The threshold value is determined by the sensitivity given by the user, where the sensitivity is the percentage of occurrences to be counted as present in the resulting presence/absence layer. For example, with a sensitivity of 0.9 or 90%, The threshold is the suitability value where 90% of the occurrence points fall on grid cells with suitability scores above that value.
The function can try multiple different sensitivity levels, and for each will output a presence/absence spatRaster stack, as well as the specificity (proportion of psuedoabsences correctly considered absent), the true skill statistic (TSS), which is a measure of how well the threshold balances commission and omission error, and the suitability value of the threshold.
The user can also downweight sensitivity when calculating TSS with the "weights" parameter. This may be important to do as 3D niche models often have a higher ratio of psuedoabsences to occurrences than 2D models.
threshold_3D(
predicted_layers,
thresholding_vals,
maxent_df,
coord_df,
weights = NULL
)A spatRaster stack of suitability layers where each layer corresponds to a depth slice
A vector of sensitivity levels for creating different thresholds and presence/absence layers. If one wanted to test sensitivities of 90%, and 95%, this would be input as c(0.9, 0.95)
'data.frame' where first column is a vector of presences named "p" containing 1's and 0's. Each row represents a cell in the spatRaster volume with an x, y, z coordinate, and 1's are presences while 0's are absences, or background points. Other columns are environmental variable values extracted at the occurrence and background points.
A dataframe containing the longitude, latitude, and depth for each cell in maxent_df, named "longitude," "latitude," and "depth"
a numeric giving what the sensitivity should be downweighted by. If no value is given, TSS will be calculated according to its original formula
a 'list' with two components:
$threshold_layers, a spatRaster stack of presence absence rasters thresholded at the sensitivity with the highest TSS
$tss_results, a 'data.frame' containing the specificity, TSS, and suitability score for each sensitivity given in thresholding_vals
Allouche O, Tsoar A, and Kadmon R. 2006. Assessing the accuracy of species distribution models: prevalence, kappa and the true skill statistic (TSS). Journal of Applied Ecology, 43, 1223-1232.
library(terra)
library(dplyr)
# creating list of spatRaster stacks where each element is a depth slice
r1_d1 <- rast(ncol = 100, nrow = 100)
set.seed(0)
values(r1_d1) <- sample(c(1:100), size = 1000, replace = TRUE)
#> Warning: [setValues] values were recycled
r2_d1 <- rast(ncol = 100, nrow = 100)
set.seed(0)
values(r2_d1) <- sample(c(1:1000), size = 1000, replace = TRUE)
#> Warning: [setValues] values were recycled
r1_d2 <- r1_d1
values(r1_d2) <- values(r1_d1)+10
r2_d2 <- r2_d1
values(r2_d2) <- values(r2_d1)+10
d1 <- c(r1_d1, r2_d1)
names(d1) <- c("valsr1", "valsr2")
d2 <- c(r1_d2, r2_d2)
names(d2) <- c("valsr1", "valsr2")
envlist <- list(d1, d2)
# creating occs and bgs
set.seed(0)
occs <- sample(c(1:nrow(crds(envlist[[1]][[1]]))), size = 50, replace = FALSE)
bgs <- sample(c(1:nrow(crds(envlist[[1]][[1]]))), size = 500, replace = FALSE)
occs_d1 <- crds(envlist[[1]][[1]])[occs[1:25],]
occs_d2 <- crds(envlist[[2]][[1]])[occs[26:50],]
bg_d1 <- crds(envlist[[1]][[1]])[bgs[1:250],]
bg_d2 <- crds(envlist[[1]][[1]])[bgs[251:500],]
# extracting at occs and bgs
occ_valsr1_d1 <- extract(envlist[[1]][[1]], occs_d1)
occ_valsr1_d2 <- extract(envlist[[2]][[1]], occs_d2)
occ_valsr2_d1 <- extract(envlist[[1]][[2]], occs_d1)
occ_valsr2_d2 <- extract(envlist[[2]][[2]], occs_d2)
occ_valsr1 <- rbind(occ_valsr1_d1, occ_valsr1_d2)
occ_valsr2 <- rbind(occ_valsr2_d1, occ_valsr2_d2)
bg_valsr1_d1 <- extract(envlist[[1]][[1]], bg_d1)
bg_valsr1_d2 <- extract(envlist[[2]][[1]], bg_d2)
bg_valsr2_d1 <- extract(envlist[[1]][[2]], bg_d1)
bg_valsr2_d2 <- extract(envlist[[2]][[2]], bg_d2)
bg_valsr1 <- rbind(bg_valsr1_d1, bg_valsr1_d2)
bg_valsr2 <- rbind(bg_valsr2_d1, bg_valsr2_d2)
valsr1 <- rbind(occ_valsr1, bg_valsr1)
valsr2 <- rbind(occ_valsr2, bg_valsr2)
p1 <- rep(1, times = 50)
p0 <- rep(0, times = 500)
p <- c(p1, p0)
maxdf <- data.frame(p, valsr1, valsr2)
# creating coord_df
coord_df <- rbind(occs_d1, occs_d2, bg_d1, bg_d2)
z <- c(rep(1, times = 25), rep(2, times = 25), rep(1, times = 250),
rep(2, times = 250))
coord_df <- cbind(coord_df, z)
colnames(coord_df) <- c("longitude", "latitude", "depth")
# creating suitability rasters
suitd1 <- d1[[1]]
values(suitd1) <- runif(min = 0, max = 1, n = length(values(suitd1)))
suitd2 <- d2[[1]]
values(suitd2) <- runif(min = 0, max = 1, n = length(values(suitd2)))
suit <- c(suitd1, suitd2)
# here's the function
result <- threshold_3D(predicted_layers = suit, thresholding_vals = c(0.9, 0.95),
maxent_df = maxdf, coord_df = coord_df, weights = 2/3)