```
### Time-specific model calibration and selection
As a final step, we will calibrate time-specific niche models using
minimum volume ellipsoids. To achieve this, we first select the
environmental variables using the function `tenm::correlation_finder`.
This function filters variables according to a correlation threshold,
which is important to avoid issues related to collinearity.
``` r
varcorrs <- tenm::correlation_finder(environmental_data =
abex$env_data[,-ncol(abex$env_data)],
method = "spearman",
threshold = 0.8,
verbose = FALSE)
#> Warning in stats::cor(environmental_data, method = method): the standard
#> deviation is zero
# Selected variables
vars2fit <- varcorrs$descriptors
print(vars2fit)
#> [1] "bio_01" "bio_02" "bio_03" "bio_04" "bio_07" "bio_12" "bio_14" "bio_15"
#> [9] "bio_17"
```
Now, we use the function `tenm::tenm_selection` to calibrate the
time-specific niche models. This function uses the background object
(here, the `abbg` object) as input. To parametrize the function, we need
to specify the omission rate criteria to be used to select the models,
the proportion of points to be included in the ellipsoid model
(*ellipsoid_level* parameter), the names of the modeling layers
(*vars2fit* parameter), a numeric vector indicating the number of
dimensions used to build ellipsoid models (*vars2fit* parameter) a
logical argument that determines whether to use the partial ROC test or
not, the random percent of data to be used for the bootstrap of the
partial ROC test (*RandomPercent* parameter), the number of iterations
of the partial ROC test (*NoOfIteration* parameter), a logical argument
to specify whether to run the calibration process in parallel and the
number of cores used in the parallel process (*parallel* parameter).
``` r
mod_sel <- tenm::tenm_selection(this_species = abbg,
omr_criteria =0.1,
ellipsoid_level=0.975,
vars2fit = vars2fit,
nvars_to_fit=c(2,3,4,5,6,7),
proc = T,
RandomPercent = 50,
NoOfIteration=1000,
parallel=TRUE,
n_cores=4)
#> -----------------------------------------------------------------------------------------
#> **** Starting model selection process ****
#> -----------------------------------------------------------------------------------------
#>
#> A total number of 36 models will be created for combinations of 9 variables taken by 2
#>
#> A total number of 84 models will be created for combinations of 9 variables taken by 3
#>
#> A total number of 126 models will be created for combinations of 9 variables taken by 4
#>
#> A total number of 126 models will be created for combinations of 9 variables taken by 5
#>
#> A total number of 84 models will be created for combinations of 9 variables taken by 6
#>
#> A total number of 36 models will be created for combinations of 9 variables taken by 7
#>
#> -----------------------------------------------------------------------------------------
#> **A total number of 492 models will be tested **
#>
#> -----------------------------------------------------------------------------------------
#> Doing calibration from model 1 to 100 in process 1
#>
#> Doing calibration from model 101 to 200 in process 2
#>
#> Doing calibration from model 201 to 300 in process 3
#>
#> Doing calibration from model 301 to 400 in process 4
#>
#> Doing calibration from model 401 to 492 in process 5
#>
#> Finishing calibration of models 1 to 100
#>
#> Finishing calibration of models 101 to 200
#>
#> Finishing calibration of models 201 to 300
#>
#> Finishing calibration of models 301 to 400
#>
#> Finishing calibration of models 401 to 492
#>
#> Finishing...
#>
#> -----------------------------------------------------------------------------------------
#> 244 models passed omr_criteria for train data
#> 27 models passed omr_criteria for test data
#> 27 models passed omr_criteria for train and test data
```
We fitted 492 models, from which 27 passed our selection criteria. Let’s
explore the `mod_sel` object.
``` r
names(mod_sel)
#> [1] "temporal_df" "sp_date_var" "lon_lat_vars" "layers_ext" "env_bg"
#> [6] "mods_table"
```
It has five slots. We can obtain the table of results of the selection
process by calling the `mods_table` slot.
``` r
head(mod_sel$mods_table,27)
#> fitted_vars nvars om_rate_train non_pred_train_ids
#> 1 bio_01,bio_02,bio_04,bio_07 4 0.06250 18,31
#> 2 bio_01,bio_02,bio_03,bio_04 4 0.06250 18,31
#> 3 bio_01,bio_03,bio_04,bio_07 4 0.06250 18,31
#> 4 bio_01,bio_04,bio_07,bio_12 4 0.09375 21,28,31
#> 5 bio_01,bio_02,bio_03,bio_07 4 0.03125 18
#> 6 bio_01,bio_04,bio_07 3 0.06250 18,31
#> 7 bio_01,bio_03,bio_04,bio_12 4 0.09375 18,21,28
#> 8 bio_01,bio_03,bio_04 3 0.06250 3,18
#> 9 bio_01,bio_04 2 0.03125 18
#> 10 bio_01,bio_02,bio_04 3 0.09375 3,18,31
#> 11 bio_01,bio_02 2 0.09375 3,18,31
#> 12 bio_01,bio_07 2 0.06250 18,31
#> 13 bio_01,bio_03 2 0.06250 3,18
#> 14 bio_01,bio_03,bio_12 3 0.06250 3,18
#> 15 bio_02,bio_04,bio_07,bio_12 4 0.06250 21,28
#> 16 bio_02,bio_03,bio_07,bio_12 4 0.06250 21,28
#> 17 bio_02,bio_03,bio_04,bio_12 4 0.06250 21,28
#> 18 bio_01,bio_07,bio_12 3 0.06250 18,31
#> 19 bio_02,bio_03,bio_04,bio_07 4 0.03125 3
#> 20 bio_04,bio_07,bio_12 3 0.06250 21,28
#> 21 bio_04,bio_07 2 0.03125 3
#> 22 bio_04,bio_12 2 0.06250 10,21
#> 23 bio_03,bio_04 2 0.03125 3
#> 24 bio_02,bio_07,bio_12 3 0.09375 3,21,28
#> 25 bio_02,bio_03,bio_12 3 0.09375 3,21,28
#> 26 bio_03,bio_07,bio_12 3 0.09375 3,21,28
#> 27 bio_07,bio_12 2 0.03125 28
#> om_rate_test non_pred_test_ids bg_prevalence pval_bin pval_proc
#> 1 0 0.4706024 0 0
#> 2 0 0.4554527 0 0
#> 3 0 0.4639894 0 0
#> 4 0 0.4627871 0 0
#> 5 0 0.4088013 0 0
#> 6 0 0.4648311 0 0
#> 7 0 0.4625466 0 0
#> 8 0 0.4655525 0 0
#> 9 0 0.4835878 0 0
#> 10 0 0.4768546 0 0
#> 11 0 0.4823855 0 0
#> 12 0 0.5001804 0 0
#> 13 0 0.4995792 0 0
#> 14 0 0.4844295 0 0
#> 15 0 0.6570879 0 0
#> 16 0 0.5706385 0 0
#> 17 0 0.6412168 0 0
#> 18 0 0.4962126 0 0
#> 19 0 0.6100757 0 0
#> 20 0 0.6847421 0 0
#> 21 0 0.7537574 0 0
#> 22 0 0.7220151 0 0
#> 23 0 0.7382470 0 0
#> 24 0 0.6783696 0 0
#> 25 0 0.6816160 0 0
#> 26 0 0.7000120 0 0
#> 27 0 0.7766021 0 0
#> env_bg_paucratio env_bg_auc mean_omr_train_test rank_by_omr_train_test
#> 1 1.520879 0.7985825 0.031250 11
#> 2 1.505433 0.7996963 0.031250 7
#> 3 1.500235 0.7974813 0.031250 8
#> 4 1.467522 0.7828762 0.046875 22
#> 5 1.465208 0.7519787 0.015625 1
#> 6 1.462961 0.7858313 0.031250 9
#> 7 1.448599 0.7549675 0.046875 21
#> 8 1.431383 0.7479400 0.031250 10
#> 9 1.429927 0.7371300 0.015625 2
#> 10 1.424194 0.7454525 0.046875 23
#> 11 1.421552 0.7490962 0.046875 24
#> 12 1.409348 0.7345675 0.031250 15
#> 13 1.409145 0.7019625 0.031250 14
#> 14 1.382147 0.6807450 0.031250 12
#> 15 1.379969 0.7115000 0.031250 18
#> 16 1.374939 0.6537675 0.031250 16
#> 17 1.373252 0.7036712 0.031250 17
#> 18 1.365017 0.7032725 0.031250 13
#> 19 1.329613 0.6943437 0.015625 3
#> 20 1.318842 0.6801412 0.031250 19
#> 21 1.307749 0.6769513 0.015625 5
#> 22 1.288655 0.6388300 0.031250 20
#> 23 1.267623 0.6472100 0.015625 4
#> 24 1.261588 0.6146950 0.046875 25
#> 25 1.253496 0.5970800 0.046875 26
#> 26 1.240056 0.5891425 0.046875 27
#> 27 1.193064 0.5648687 0.015625 6
#> rank_omr_aucratio
#> 1 1
#> 2 2
#> 3 3
#> 4 4
#> 5 5
#> 6 6
#> 7 7
#> 8 8
#> 9 9
#> 10 10
#> 11 11
#> 12 12
#> 13 13
#> 14 14
#> 15 15
#> 16 16
#> 17 17
#> 18 18
#> 19 19
#> 20 20
#> 21 21
#> 22 22
#> 23 23
#> 24 24
#> 25 25
#> 26 26
#> 27 27
```
### Projecting time-specific niche models
To project the models, we use the `predict` method. Here, we will
project one of the selected models using the environmental layers of
2016. Also, we project it using layers from a period that comprehends
1970-2000.
``` r
env_layers_2016 <- list.dirs(tempora_layers_dir,
recursive = FALSE)[32]
suit_2016 <- predict(mod_sel,
model_variables = c("bio_01","bio_03","bio_12"),
layers_path =env_layers_2016 ,
layers_ext = ".tif$")
#> | | | 0% | |======================================================================| 100%
```
Now for the period that comprehends 1970-2000.
``` r
layers_70_00_dir <- system.file("extdata/bio_1970_2000",package = "tenm")
suit_1970_2000 <- predict(mod_sel,
model_variables = c("bio_01","bio_03","bio_12"),
layers_path = layers_70_00_dir,
layers_ext = ".tif$")
#> | | | 0% | |======================================================================| 100%
```
Lets see the predictions in geographic space
``` r
par(mfrow=c(1,2), mar=c(4,4,2,2))
terra::plot(suit_2016, main="Prediction for 2016")
terra::plot(suit_1970_2000, main="Prediction for 1970-2000")
```
### Comparing time-specific niche model vs. standard niche model
The following lines of code show the differences of a time-specific
niche model and a standard niche model.
``` r
layers_70_00_dir <- system.file("extdata/bio_1970_2000",package = "tenm")
layers_70_00_path <- list.files(layers_70_00_dir,
pattern = ".tif$",full.names = TRUE)
# Extract environmental information
elayers_70_00 <- terra::rast(layers_70_00_path)
e_trad <- terra::extract(elayers_70_00,
ab_1[,c("decimalLongitude","decimalLatitude")])
rgl::view3d(theta = 0, phi = -60,fov=120, zoom = 0.7)
tenm::plot_ellipsoid(x = e_trad$bio_01,y=e_trad$bio_03,z=e_trad$bio_12,
col = "#1B9E77",
xlab = "Bio 1",
ylab = "Bio 3",
zlab = "Bio 12",)
tenm::plot_ellipsoid(x = abbg$temporal_df$bio_01,
y = abbg$temporal_df$bio_03,
z = abbg$temporal_df$bio_12,
col = "#E7298A",
add = TRUE)
```
Note that both ellipsoids differ in size and shape. In standard approach
(green ellipsoid), we can see an sub-estimation of the environmental
values where the intrinsic growth rate might be positive.
