library(knitr) knitr::opts_chunk$set(echo = TRUE, warning = FALSE, message = FALSE, cache = FALSE, fig.width = 6, fig.height = 6, fig.show = "hold")
Start with the necessary packages and seed for the vignette.
loadedPackages <- c("envi", "RStoolbox", "sf", "spatstat.data", "spatstat.geom", "spatstat.random", "terra") invisible(lapply(loadedPackages, library, character.only = TRUE)) set.seed(1234) # for reproducibility
We use the gorillas
data and the accompanying covariate data in gorillas.extra
from the spatstat.data
package on CRAN. These data are locations of nesting sites of gorillas in the Kagwene Gorilla Sanctuary in Cameroon. A detailed description and analysis of the data are reported in Funwi-Gabga and Mateu (2012). The authors used a kernel density-based smoothing technique to detect hot spots of nesting in the park. Here, we use another kernel density-based smoothing technique to detect hot spots of nesting within the covariate information (i.e., the gorilla ecological niche) and then predict where these hot spots are located within the park.
We start by importing the two covariate data of class im
:
slopeangle <- spatstat.data::gorillas.extra$slopeangle waterdist <- spatstat.data::gorillas.extra$waterdist
Center and scale the covariate data.
slopeangle$v <- scale(slopeangle) waterdist$v <- scale(waterdist)
Convert the covariate data to class SpatRaster
.
slopeangle_raster <- terra::rast(slopeangle) waterdist_raster <- terra::rast(waterdist)
Add appropriate marks to the gorillas
data from spatstat.data
package. These points are considered our "presence" locations.
presence <- spatstat.geom::unmark(spatstat.data::gorillas) spatstat.geom::marks(presence) <- data.frame("presence" = rep(1, presence$n), "lon" = presence$x, "lat" = presence$y) spatstat.geom::marks(presence)$slopeangle <- slopeangle[presence] spatstat.geom::marks(presence)$waterdist <- waterdist[presence]
Randomly draw points from the study area and add the appropriate marks. These points are considered our "(pseudo-)absence" locations.
absence <- spatstat.random::rpoispp(0.00004, win = slopeangle) spatstat.geom::marks(absence) <- data.frame("presence" = rep(0, absence$n), "lon" = absence$x, "lat" = absence$y) spatstat.geom::marks(absence)$slopeangle <- slopeangle[absence] spatstat.geom::marks(absence)$waterdist <- waterdist[absence]
rm(waterdist, slopeangle) # conserve memory
Combine the presence (n = r presence$n
) and absence (r absence$n
) locations into one object of class data.frame
and reorder the features required for the lrren
function in the envi
package:
obs_locs <- spatstat.geom::superimpose(absence, presence, check = FALSE) spatstat.geom::marks(obs_locs)$presence <- as.factor(spatstat.geom::marks(obs_locs)$presence) spatstat.geom::plot.ppp(obs_locs, which.marks = "presence", main = "Gorilla nesting sites (red-colored)\nPseudo-absence locations (blue-colored)", cols = c("#0000CD","#8B3A3A"), pch = 1, axes = TRUE, ann = TRUE) obs_locs <- spatstat.geom::marks(obs_locs) obs_locs$id <- seq(1, nrow(obs_locs), 1) obs_locs <- obs_locs[ , c(6, 2, 3, 1, 4, 5)]
Extract the prediction locations within the study area from one of the covariates.
predict_xy <- terra::crds(slopeangle_raster) predict_locs <- as.data.frame(predict_xy) predict_locs$slopeangle <- terra::extract(slopeangle_raster, predict_xy)[ , 1] predict_locs$waterdist <- terra::extract(waterdist_raster, predict_xy)[ , 1]
rm(waterdist_raster, slopeangle_raster) # conserve memory
Run the lrren
function within the envi
package. We use the default settings, except we want to predict the ecological niche within the study area (predict = TRUE
), so we conduct k-fold cross-validation model fit diagnostics (cv = TRUE
) by undersampling absence locations to balance the prevalence (0.5) within all testing data sets (balance = TRUE
).
start_time <- Sys.time() # record start time test_lrren <- envi::lrren(obs_locs = obs_locs, predict_locs = predict_locs, predict = TRUE, cv = TRUE, balance = TRUE) end_time <- Sys.time() # record end time lrren_time <- end_time - start_time # calculate duration of lrren() example
A single run of the lrren
function took approximately r round(lrren_time, digits = 1)
seconds on a machine with the features listed at the end of the vignette.
We display the estimated ecological niche within a space of Covariate 1 by Covariate 2 using the plot_obs
function. We use the default two-tailed alpha level (alpha = 0.05
) and the default colors where the yellow color denotes areas with covariate data combinations where we have sparse observations. As expected, extreme values of the log relative risk surface are located near the edges of the surface; however, these areas are highly variable and are not statistically significant based on an asymptotic normal assumption. The default color key for the log relative risk surface hides the heterogeneity closer to the null value (zero). Therefore, we limit the color key for the log relative risk surface to (-1, 1).
envi::plot_obs(test_lrren, lower_lrr = -1, upper_lrr = 1)
We observe two areas of positive log relative risk, and the center of these areas are statistically significant, suggesting two ecological niches of gorilla nesting sites compared to pseudo-absence points drawn randomly from within the park (based on only two covariates, slope gradient and distance from water). One niche is a combination of flat to moderate (about 10 - 30 degrees) slope gradient and moderate to far (about 200 - 400 meters) distance from water, which may correspond to the top of ridges. The second niche is a combination of moderate (30 - 40 degrees) slope gradient and moderate (about 100 - 200 meters) distance from water, which may correspond to within valleys.
We display the estimated ecological niche predicted to the study area within geographic space using the plot_predict
function. We use the default two-tailed alpha level (alpha = 0.05
) and the default colors where the yellow color denotes areas with covariate data combinations where we have sparse observations. We limit the color key for the log relative risk prediction to (-1, 1).
envi::plot_predict(test_lrren, cref0 = "EPSG:32632", cref1 = "EPSG:4326", lower_lrr = -1, upper_lrr = 1)
The two estimated ecological niches are located in many small regions throughout the park, reflected by the large spatial heterogeneity in the two covariates. For example, the tops of ridges are the farthest distance from water and are located in many areas throughout the park. Importantly, gorilla nesting sites were not observed in many of these areas, but this prediction represents areas with combinations of covariates that match (or are similar) to the combinations occurring at observed locations.
This is an example of a scale mismatch. The scale of the gorilla nesting site presence is a broad, elliptical-shaped area in the northwest region of the park. Spatial interpolation of the gorilla nesting sites (i.e., not considering covariate information) can be seen in Figure 5a within the original study by Funwi-Gabga and Mateu (2012). The two covariates (slope gradient and distance from water) vary considerably within this area, and we observe a full range of covariate values within the larger gorilla presence area. Our approach is a version of environmental interpolation and considers covariate information when predicting the spatial distribution of gorilla nesting sites. When the gorilla nesting sites are arranged in covariate space, the smoothing bandwidth of our kernel density-based approach creates a smoother relative risk surface than the covariate surfaces themselves because the gorilla nesting sites are clustering in particular combinations of the two covariates (i.e., ecological niche) that are scattered throughout the park in geographic space.
We display the internal 10-fold cross-validation diagnostics using the plot_cv
function. We use the default two-tailed alpha level (alpha = 0.05
) and our prevalence is fairly balanced at r presence$n/(presence$n + absence$n)
.
envi::plot_cv(test_lrren)
rm(absence, presence, test_lrren) # conserve memory
The log relative risk estimate accurately predicts about 60% of the gorilla nesting sites, which is better than chance (50%) but not a large improvement. The pseudo-absence locations are drawn at random throughout the park and are located in areas estimated with higher log relative risk values, which reduces the prediction performance of the model. Using observed absences instead of pseudo-absences may improve the prediction of the spatial distribution of gorilla nesting sites.
The choice in covariates is critical for ecological niche models and especially for lrren
because we are limited to two covariates. The goal, if possible, is to discover covariates that separate presence locations (i.e., nesting sites) from absence locations. Because our pseudo-absence locations are randomly drawn within the park it will be challenging to completely separate presence and absence locations in covariate space. One approach to include more than two covariates is to conduct an orthogonal linear transformation on multiple covariates, such as a Principal Component Analysis (PCA). Here, we can use the first two components of a PCA of all seven available covariates in the gorillas.extra
data in the spatstat.data
package, which include:
We can use the rasterPCA
function within the RStoolbox
package to conduct a PCA of multiple spatial data layers. We start by centering and scaling the numeric-valued layers. We also center the factor-valued layer 'aspect' at "North" and we group flatter slope types together in the factor-valued layer 'slope type.' NOTE: Using categorical (discrete) variables in PCA is requires more consideration than demonstrated in this vignette (see: Kolenikov & Angeles (2009) for more details).
aspect <- spatstat.data::gorillas.extra$aspect # class factor elevation <- spatstat.data::gorillas.extra$elevation # class integer heat <- spatstat.data::gorillas.extra$heat # class factor slopeangle <- spatstat.data::gorillas.extra$slopeangle # class numeric slopetype <- spatstat.data::gorillas.extra$slopetype # class factor vegetation <- spatstat.data::gorillas.extra$vegetation # class factor waterdist <- spatstat.data::gorillas.extra$waterdist # class numeric # Center and scale numeric elevation$v <- scale(elevation$v) slopeangle$v <- scale(slopeangle$v) waterdist$v <- scale(waterdist$v) # Create rasters aspect <- terra::rast(aspect) elevation <- terra::rast(elevation) heat <- terra::rast(heat) slopeangle <- terra::rast(slopeangle) slopetype <- terra::rast(slopetype) vegetation <- terra::rast(vegetation) waterdist <- terra::rast(waterdist) # Reorder aspect to center by "N" instead of "S" terra::values(aspect) <- factor(terra::values(aspect), levels = c("5", "6", "7", "8", "1", "2", "3", "4")) # Reorder slope types to order flatter types next to each other terra::values(slopetype) <- factor(terra::values(slopetype), levels = c("1", "6", "3", "2", "4", "5")) # Stack of SpatRasters park <- c(aspect, elevation, heat, slopeangle, slopetype, vegetation, waterdist) # Principal Component Analysis pca_park <- RStoolbox::rasterPCA(park) summary(pca_park$model) # PCA components
rm(aspect, elevation, heat, slopeangle, slopetype, vegetation, waterdist, park) # conserve memory
The first two components of the PCA explain almost 64% of the variation across the seven covariates.
pca_park$model$loadings # PCA loadings
The loadings of the first component is almost entirely the centered 'aspect' variable. The loadings of the second component is a combination of centered 'elevation,' 'vegetation type,' centered 'distance from water,' and centered 'slope type.'
# Extract Bands from PCA pca_bands <- pca_park$map pc1 <- pca_bands[[1]] # PC1 pc2 <- pca_bands[[2]] # PC2 pc1 <- envi:::as.im.SpatRaster(terra::rast(pc1)) # convert to class 'im' pc2 <- envi:::as.im.SpatRaster(terra::rast(pc2)) # convert to class 'im' spatstat.geom::plot.im(pc1, main = 'Principal Component 1\nprimarily aspect (centered at "North")', ann = TRUE, axes = TRUE) spatstat.geom::plot.im(pc2, main = 'Principal Component 2\ncombination of elevation, vegetation type,\ndistance from water, and slope type', ann = TRUE, axes = TRUE)
rm(pca_park) # conserve memory
For the first component (centered aspect), we can observe northern (i.e., North, North East, North West) aspects are primarily located in the northern section of the park. For the second component (combination), we can see a large cluster of middle loading values (between 0 and 2) in the center of the park surrounded by lower loading values (< 0) and a few smaller clusters of high loading values (> 2) and low loading values (< -2) within the large middle loading value cluster.
We prepare inputs for a new lrren
run similar to the first example above. For consistency, we also reset the random number generator to select similar pseudo-absence locations as the example above.
set.seed(1234) # for similar locations as above example presence <- spatstat.geom::unmark(spatstat.data::gorillas) spatstat.geom::marks(presence) <- data.frame("presence" = rep(1, presence$n), "lon" = presence$x, "lat" = presence$y) spatstat.geom::marks(presence)$pc1 <- pc1[presence] spatstat.geom::marks(presence)$pc2 <- pc2[presence] absence <- spatstat.random::rpoispp(0.00004, win = spatstat.data::gorillas.extra$aspect) spatstat.geom::marks(absence) <- data.frame("presence" = rep(0, absence$n), "lon" = absence$x, "lat" = absence$y) spatstat.geom::marks(absence)$pc1 <- pc1[absence] spatstat.geom::marks(absence)$pc2 <- pc2[absence] obs_locs <- spatstat.geom::superimpose(absence, presence, check = FALSE) spatstat.geom::marks(obs_locs)$presence <- as.factor(spatstat.geom::marks(obs_locs)$presence) obs_locs <- spatstat.geom::marks(obs_locs) obs_locs$id <- seq(1, nrow(obs_locs), 1) obs_locs <- obs_locs[ , c(6, 2, 3, 1, 4, 5)] predict_xy <- terra::crds(terra::rast(pca_bands[[1]])) predict_locs <- as.data.frame(predict_xy) predict_locs$pc1 <- terra::extract(terra::rast(pca_bands[[1]]), predict_xy)[ , 1] predict_locs$pc2 <- terra::extract(terra::rast(pca_bands[[2]]), predict_xy)[ , 1]
rm(pca_bands) # conserve memory
With the two principal components, we re-run the lrren
function.
start_time <- Sys.time() # record start time pca_lrren <- envi::lrren(obs_locs = obs_locs, predict_locs = predict_locs, predict = TRUE, cv = TRUE, balance = TRUE) end_time <- Sys.time() # record end time pca_time <- end_time - start_time # calculate duration of lrren() example
A single run of the lrren
function took approximately r round(pca_time, digits = 1)
seconds on a machine with the features listed at the end of the vignette.
We display the estimated ecological niche within a space of Principal Component 1 by Principal Component 2 using the plot_obs
function. We use the default two-tailed alpha level (alpha = 0.05
) and the default colors where the yellow color denotes areas with covariate data combinations where we have sparse observations. We limit the color key for the log relative risk surface to (-1, 1).
envi::plot_obs(pca_lrren, lower_lrr = -1, upper_lrr = 1)
We display the estimated ecological niche predicted to the study area within geographic space using the plot_predict
function. We use the default two-tailed alpha level (alpha = 0.05
), and we limit the color key for the log relative risk prediction to (-1, 1).
envi::plot_predict(pca_lrren, cref0 = "EPSG:32632", cref1 = "EPSG:4326", lower_lrr = -1, upper_lrr = 1)
We display the internal 10-fold cross-validation diagnostics using the plot_cv
function. We use the default two-tailed alpha level (alpha = 0.05
).
envi::plot_cv(pca_lrren)
rm(pca_lrren, obs_locs, predict_locs) # conserve memory
Based on only the first two components of the Principal Component Analysis of seven covariates, we detected one ecological niche of gorillas compared to pseudo-absence points drawn randomly from within the park. Presence and absences appear separated more by Principal Component 2 (loading values between 1 and 4) than Principal Component 1. This is reflected in geographic space where the gorilla niche is located in an area similar to the large cluster of middle loading values of Principal Component 2, which is located in the central, northwestern section of the park. The log relative risk estimate accurately predicted about 75% of the gorilla nesting sites, which is markedly improved from our first example above (60%).
The two components captured about 64% of the variation across the seven covariates, which may not be important in separating presence locations from absence locations within covariate space. Further studies can assess different combinations of principal components (e.g., Principal Component 2 and Principal Component 3) or include another yet unavailable covariate that may be ecologically important for gorillas.
The original study by Funwi-Gabga and Mateu (2012) used kernel density estimation techniques in geographic space (see: Figure 5a in reference). Our method applies a kernel density estimation technique in covariate space and shows results similar to the inhomogeneous spatial point process model results of the original study (see: Figure 10a in reference), which incorporates the same covariate information to predict the spatial distribution of gorilla nesting sites. Our approach predicted more nesting sites in the western section of the park than in the original study.
Now we perform a sensitivity analysis of our ecological niche model. For example, let us assume the investigators observed nesting sites at a distance, and there is some spatial uncertainty in the exact coordinates of the nesting sites. Let us also assume there is more uncertainty about where the nesting sites of the 'major' gorilla group were located (within 100 meters) than the 'minor' gorilla group (within 50 meters). We can examine the influence of this type of uncertainty on our estimated ecological niche model and its predicted spatial distribution in the park. We assume there is no spatial uncertainty in the spatial location of pseudo-absence points (0.1 meters). The three groups must be categorized within a new feature named levels
within the observation data.
We start by preparing the observation data. Here, the data must be a marked planar point pattern of class 'ppp'.
# Preserve 'group' feature from 'gorillas' data within {spatstat.data} ## Assign as a new mark named 'level' spatstat.geom::marks(presence)$levels <- spatstat.geom::marks(spatstat.data::gorillas)$group # Assign a third 'level' for the pseudo-absences spatstat.geom::marks(absence)$levels <- "none" # Combine obs_locs <- spatstat.geom::superimpose(absence, presence, check = FALSE) # Set variables as class 'factor' spatstat.geom::marks(obs_locs)$presence <- as.factor(spatstat.geom::marks(obs_locs)$presence) spatstat.geom::marks(obs_locs)$levels <- as.factor(spatstat.geom::marks(obs_locs)$levels) # Create 'id' feature spatstat.geom::marks(obs_locs)$id <- seq(1, obs_locs$n, 1) # Reorder and drop the two covariate features spatstat.geom::marks(obs_locs) <- spatstat.geom::marks(obs_locs)[ , c(7, 2, 3, 1, 6)]
rm(absence, presence) # conserve memory
The two covariate values will be assigned to all points in every iteration and are drawn from a list of objects of class 'im'. Here, we use the same variables as above (slope gradient and distance from water), centered and scaled.
ims <- spatstat.geom::solist(pc1,pc2) names(ims) <- c("Component_1", "Component_2")
rm(pc1, pc2) # conserve memory
Run the perlrren
function within the envi
package. We use the default settings and 10 simulated iterations. The radii
argument is a string of numeric values in the order of the levels
feature of class 'factor' within the obs_locs
object.
n_sim <- 100 start_time <- Sys.time() # record start time test_perlrren <- envi::perlrren(obs_ppp = obs_locs, covariates = ims, radii = c(100, 50, 0.1), n_sim = n_sim) end_time <- Sys.time() # record end time perlrren_time <- end_time - start_time # calculate duration of perlrren() example
rm(obs_locs, ims) # conserve memory
A single (non-parallel) run of the perlrren
function with r n_sim
iterations took approximately r round(perlrren_time, digits = 1)
minutes on a machine with the features listed at the end of the vignette. In practice, a larger number of simulated iterations would be recommended (e.g., n = 10,000).
We display the summary statistics from the sensitivity analysis with the plot_perturb
function. Here, we use the default settings, except we limit the color key for the log relative risk surface and log relative risk prediction to (-1, 1). We only display the plots pertaining to the proportion of iterations with a significant p-value, both in covariate space and predicted into geographic space (i.e., the park).
envi::plot_perturb(test_perlrren, cov_labs = c("Component 1", "Component 2"), cref0 = "EPSG:32632", cref1 = "EPSG:4326", lower_lrr = -1, upper_lrr = 1)
rm(test_perlrren) # conserve memory
The uncertainty in the spatial coordinates of gorilla nesting sites greatly influences the estimated ecological niche of gorillas (based on the two principal components) compared to pseudo-absence points drawn randomly from within the park. The significant hot spot (i.e., ecological niche) within the covariate space diminished in size when adding this type of spatial uncertainty to the model. The geographic location of areas the model cannot distinguish as suitable or unsuitable for gorilla nesting locations is in the center of the park or the southern edge of the secondary forest with a southern facing aspect.
sessionInfo()
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.