Start with the necessary packages and seed for the vignette.
Data preparation for lrren
We start by importing the two covariate data of class
im:
- Slope gradient (degrees)
- Distance from water (meters)
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 RasterLayer.
slopeangle_raster <- raster::raster(slopeangle)
waterdist_raster <- raster::raster(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]
Combine the presence (n = 647) and absence (769) locations into one
object of class data.frame and reorder the features
required for the lrren function in the envi
package:
- ID
- X-coordinate
- Y-coordinate
- Presence (binary)
- Covariate 1
- Covariate 2
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_locs <- data.frame(raster::rasterToPoints(slopeangle_raster))
predict_locs$layer2 <- raster::extract(waterdist_raster, predict_locs[, 1:2])
lrren with two covariates
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 <- 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
2.8 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).
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).
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 0.4569209.
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.
lrren with two principal components
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:
- Aspect
- Elevation (meters)
- Heat index
- Slope gradient
- Slope position
- Vegetation type
- Distance from water
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 <- raster::raster(aspect)
elevation <- raster::raster(elevation)
heat <- raster::raster(heat)
slopeangle <- raster::raster(slopeangle)
slopetype <- raster::raster(slopetype)
vegetation <- raster::raster(vegetation)
waterdist <- raster::raster(waterdist)
# Reorder aspect to center by "N" instead of "S"
raster::values(aspect) <- factor(values(aspect),
levels = c("5", "6", "7", "8", "1", "2", "3", "4"))
# Reorder slope types to order flatter types next to each other
raster::values(slopetype) <- factor(values(slopetype),
levels = c("1", "6", "3", "2", "4", "5"))
# Save names for RasterLayers
names(aspect) <- "aspect"
names(elevation) <- "elevation"
names(heat) <- "heat"
names(slopeangle) <- "slopeangle"
names(slopetype) <- "slopetype"
names(vegetation) <- "vegetation"
names(waterdist) <- "waterdist"
# RasterStack
park <- raster::stack(aspect, elevation, heat, slopeangle, slopetype, vegetation, waterdist)
# Principal Component Analysis
pca_park <- RStoolbox::rasterPCA(park)
summary(pca_park$model) # PCA components
## Importance of components:
## Comp.1 Comp.2 Comp.3 Comp.4 Comp.5
## Standard deviation 2.1697414 1.5626701 1.5429105 1.00679778 0.97729569
## Proportion of Variance 0.3786484 0.1964060 0.1914704 0.08152761 0.07681962
## Cumulative Proportion 0.3786484 0.5750545 0.7665249 0.84805249 0.92487211
## Comp.6 Comp.7
## Standard deviation 0.81117157 0.52542744
## Proportion of Variance 0.05292315 0.02220474
## Cumulative Proportion 0.97779526 1.00000000
The first two components of the PCA explain over 57% of the variation
across the seven covariates.
pca_park$model$loadings # PCA loadings
##
## Loadings:
## Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6 Comp.7
## aspect 0.994
## elevation -0.373 -0.179 0.276 0.864
## heat 0.998
## slopeangle 0.611 0.767 -0.147
## slopetype -0.482 0.872
## vegetation -0.773 -0.436 -0.294 0.118 -0.322
## waterdist -0.155 0.678 -0.618 -0.357
##
## Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6 Comp.7
## SS loadings 1.000 1.000 1.000 1.000 1.000 1.000 1.000
## Proportion Var 0.143 0.143 0.143 0.143 0.143 0.143 0.143
## Cumulative Var 0.143 0.286 0.429 0.571 0.714 0.857 1.000
The loadings of the first component is almost entirely the centered
‘aspect’ variable. The loadings of the second component is a combination
of ‘vegetation type,’ centered ‘slope type,’ centered ‘elevation,’ and
centered ‘distance from water.’
# Extract Bands from PCA
pca_bands <- pca_park$map
pc1 <- pca_bands[[1]] # PC1
pc2 <- pca_bands[[2]] # PC2
pc1 <- spatstat.geom::as.im(pc1) # convert to class 'im'
pc2 <- spatstat.geom::as.im(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 vegetation type, slope type,\nelevation, and distance from water',
ann = TRUE,
axes = TRUE)
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 -2 and 0) in the center of the park surrounded by high loading
values (> 0) and a few smaller clusters of 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_locs <- data.frame(raster::rasterToPoints(pca_bands[[1]]))
predict_locs$layer2 <- raster::extract(pca_bands[[2]], predict_locs[, 1:2])
With the two principal components, we re-run the lrren
function.
start_time <- Sys.time() # record start time
pca_lrren <- 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 4
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).
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).
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).
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 -4 and -1) 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 73% of the gorilla nesting sites, which is markedly improved from
our first example above (60%).
The two components captured about 57% 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.
Sensitivity Analysis
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)]
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")
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 <- 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
A single (non-parallel) run of the perlrren function
with 100 iterations took approximately 1.4 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).
plot_perturb(test_perlrren,
cov_labs = c("Component 1", "Component 2"),
cref0 = "EPSG:32632",
cref1 = "EPSG:4326",
lower_lrr = -1,
upper_lrr = 1)
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.
## R version 4.2.1 (2022-06-23)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Catalina 10.15.7
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] spatstat.random_2.2-0 spatstat.geom_2.4-0 spatstat.data_2.2-0
## [4] RStoolbox_0.3.0 raster_3.5-21 maptools_1.1-4
## [7] sp_1.5-0 envi_0.1.14 knitr_1.39
##
## loaded via a namespace (and not attached):
## [1] spam_2.9-1 plyr_1.8.7 splines_4.2.1
## [4] listenv_0.8.0 ggplot2_3.3.6 digest_0.6.29
## [7] foreach_1.5.2 htmltools_0.5.3 viridis_0.6.2
## [10] fansi_1.0.3 magrittr_2.0.3 tensor_1.5
## [13] doParallel_1.0.17 ROCR_1.0-11 recipes_1.0.1
## [16] globals_0.16.0 gower_1.0.0 doFuture_0.12.2
## [19] exactextractr_0.8.1 hardhat_1.2.0 spatstat.sparse_2.1-1
## [22] colorspace_2.0-3 xfun_0.31 dplyr_1.0.9
## [25] rgdal_1.5-30 tcltk_4.2.1 jsonlite_1.8.0
## [28] spatstat_2.3-4 survival_3.3-1 iterators_1.0.14
## [31] glue_1.6.2 polyclip_1.10-0 gtable_0.3.0
## [34] ipred_0.9-13 V8_4.2.1 future.apply_1.9.0
## [37] maps_3.4.0 abind_1.4-5 scales_1.2.1
## [40] DBI_1.1.3 rngtools_1.5.2 Rcpp_1.0.9
## [43] viridisLite_0.4.1 units_0.8-0 sparr_2.2-16
## [46] spatstat.core_2.4-4 foreign_0.8-82 proxy_0.4-27
## [49] dotCall64_1.0-1 stats4_4.2.1 lava_1.6.10
## [52] prodlim_2019.11.13 pkgconfig_2.0.3 XML_3.99-0.10
## [55] nnet_7.3-17 sass_0.4.2 deldir_1.0-6
## [58] utf8_1.2.2 caret_6.0-92 tidyselect_1.1.2
## [61] rlang_1.0.4 reshape2_1.4.4 munsell_0.5.0
## [64] tools_4.2.1 cachem_1.0.6 cli_3.3.0
## [67] generics_0.1.3 pls_2.8-1 evaluate_0.15
## [70] stringr_1.4.0 fastmap_1.1.0 yaml_2.3.5
## [73] goftest_1.2-3 ModelMetrics_1.2.2.2 purrr_0.3.4
## [76] future_1.27.0 nlme_3.1-158 doRNG_1.8.2
## [79] concaveman_1.1.0 compiler_4.2.1 rstudioapi_0.13
## [82] curl_4.3.2 e1071_1.7-11 spatstat.utils_2.3-1
## [85] spatstat.linnet_2.3-2 tibble_3.1.8 bslib_0.4.0
## [88] stringi_1.7.8 highr_0.9 fields_14.1
## [91] rgeos_0.5-9 lattice_0.20-45 Matrix_1.4-1
## [94] classInt_0.4-7 vctrs_0.4.1 pillar_1.8.1
## [97] lifecycle_1.0.1 jquerylib_0.1.4 data.table_1.14.2
## [100] R6_2.5.1 KernSmooth_2.23-20 gridExtra_2.3
## [103] parallelly_1.32.1 codetools_0.2-18 MASS_7.3-57
## [106] cvAUC_1.1.4 withr_2.5.0 mgcv_1.8-40
## [109] parallel_4.2.1 terra_1.5-21 grid_4.2.1
## [112] rpart_4.1.16 timeDate_4021.104 class_7.3-20
## [115] rmarkdown_2.14 misc3d_0.9-1 sf_1.0-7
## [118] pROC_1.18.0 lubridate_1.8.0