Projecting models
Andrea Sánchez Tapia
2019-11-12
Here we test the function setup_sdmdata() and its different parametrization options. This function prepares the data to run the models with do_any() or do_many() functions. It creates an object in the workspace to be used in do_any() or do_many() and writes two csv files with metadata and sdmdata.
Loading required packages.
library(rJava)
library(raster)
#library(modleR)
library(dplyr)
#eu estou usando uma cópia local para desenvolvimento
devtools::load_all("../../1_modleR")
library(maps)
library(maptools)Selecting one species as an example
We use a standard dataset inside the package modleR. First, from example_occs object we select only data from one species Abarema langsdorffii and create one training set (70% of the data) and one test set (30% of the data) for the data.
## Creating an object with species names
especies <- names(example_occs)[1]
# Selecting only coordinates for the first species
coord1sp <- example_occs[[1]]
head(coord1sp)## sp lon lat
## 343 Abarema_langsdorffii -40.615 -19.921
## 344 Abarema_langsdorffii -40.729 -20.016
## 345 Abarema_langsdorffii -41.174 -20.303
## 346 Abarema_langsdorffii -41.740 -20.493
## 347 Abarema_langsdorffii -42.482 -20.701
## 348 Abarema_langsdorffii -40.855 -17.082dim(coord1sp)## [1] 104 3# Subsetting data into training and test
ceiling(0.7 * nrow(coord1sp))## [1] 73# Making a sample of 70% of species' records
set <- sample(1:nrow(coord1sp), size = ceiling(0.7 * nrow(coord1sp)))
# Creating training data set (70% of species' records)
train_set <- coord1sp[set,]
# Creating test data set (other 30%)
test_set <- coord1sp[setdiff(1:nrow(coord1sp),set),]Now lets the check our data points. We plot the traning and test data sets with the first axis of the environmental PCA data from the object example_vars.
# selecting only the first PCA axis
predictor <- example_vars[[1]]
# transforming the data frame with the coordinates in a spatial object
pts <- SpatialPoints(coord1sp[,c(2,3)])
# ploting environmental layer
plot(predictor, legend = F)
points(lat ~ lon, data=coord1sp)
Fitting and projection datasets
We want to fit everything in a small area around the occurrence points (generated through a “median distance buffer”) and project it to other sample rasterStacks
fit_data <- stack("../data/env/cropped_proj.tif")
proj_folder <- "../data/env/proj/proj1/"
proj_data <- list.files(proj_folder, full.names = T) %>%
stack() %>%
.[[1]]
pts <- SpatialPoints(coord1sp[,c(2,3)])
plot(!is.na(example_vars[[1]]), legend = F, add = F)
plot(fit_data[[1]], legend = F, add = T)
plot(proj_data[[1]], legend = F, add = T)
points(train_set[,2:3], col = "red", pch = 19)
points(test_set[,2:3], col = "blue", pch = 19)
Model projection
A simple model with no buffer and no projection:
# generating sdmdata
sdm_no_proj <- setup_sdmdata(species_name = especies[1],
occurrences = coord1sp[,-1],
predictors=fit_data,
models_dir="projection/fit",
clean_dupl = T,
clean_na = T)
# running models
no_proj <- do_many(species_name = especies[1],
predictors = fit_data,
models_dir = "projection/fit",
bioclim = T,
svmk = T)
no_proj.bc <- raster("projection/fit/Abarema_langsdorffii/present/partitions/bioclim_cont_Abarema_langsdorffii_1_1.tif")
no_proj.svmk <- raster("projection/fit/Abarema_langsdorffii/present/partitions/svmk_cont_Abarema_langsdorffii_1_1.tif")
plot(no_proj.bc)
maps::map(, , add = T)
plot(no_proj.svmk)
maps::map(, , add = T)
We project to a series of projections in the specified folder: "./env/proj"
## first renaming layer names from predictors to match proj files
names(fit_data) <- paste0("proj1_", 1:6)
sdm_yes_proj <- setup_sdmdata(species_name = especies[1],
occurrences = coord1sp[,-1],
predictors = fit_data,
models_dir = "projection/projtest",
clean_dupl = T,
clean_na = T,
buffer_type = "median",
partition_type = "bootstrap",
boot_n = 3)
yes_proj1 <- do_any(species_name = especies[1],
predictors = fit_data,
models_dir = "projection/projtest",
algo = "bioclim",
project_model = T,
proj_data_folder = "data/env/proj/",
write_png = T)
yes_proj <- do_many(species_name = especies[1],
predictors = fit_data,
models_dir = "projection/projtest",
bioclim = T,
domain = T,
glm = T,
mahal = F,
maxent = T,
rf = T,
svmk = T,
svme = T,
brt = T,
project_model = T,
proj_data_folder = "../data/env/proj/",
write_png = T)## Warning in randomForest.default(x, y, mtry = mtryStart, ntree = ntreeTry, :
## The response has five or fewer unique values. Are you sure you want to do
## regression?## Warning in randomForest.default(x, y, mtry = mtryCur, ntree = ntreeTry, :
## The response has five or fewer unique values. Are you sure you want to do
## regression?## -0.04878859 0.05## Warning in randomForest.default(x, y, mtry = mtryCur, ntree = ntreeTry, :
## The response has five or fewer unique values. Are you sure you want to do
## regression?## -0.1160929 0.05## Warning in randomForest.default(x, y, mtry = res[which.min(res[, 2]), 1], :
## The response has five or fewer unique values. Are you sure you want to do
## regression?## Warning in randomForest.default(x, y, mtry = mtryStart, ntree = ntreeTry, :
## The response has five or fewer unique values. Are you sure you want to do
## regression?## Warning in randomForest.default(x, y, mtry = mtryCur, ntree = ntreeTry, :
## The response has five or fewer unique values. Are you sure you want to do
## regression?## 0.02297148 0.05## Warning in randomForest.default(x, y, mtry = mtryCur, ntree = ntreeTry, :
## The response has five or fewer unique values. Are you sure you want to do
## regression?## -0.04537304 0.05## Warning in randomForest.default(x, y, mtry = res[which.min(res[, 2]), 1], :
## The response has five or fewer unique values. Are you sure you want to do
## regression?## Warning in randomForest.default(x, y, mtry = mtryStart, ntree = ntreeTry, :
## The response has five or fewer unique values. Are you sure you want to do
## regression?## Warning in randomForest.default(x, y, mtry = mtryCur, ntree = ntreeTry, :
## The response has five or fewer unique values. Are you sure you want to do
## regression?## -0.09380388 0.05## Warning in randomForest.default(x, y, mtry = mtryCur, ntree = ntreeTry, :
## The response has five or fewer unique values. Are you sure you want to do
## regression?## -0.03368778 0.05## Warning in randomForest.default(x, y, mtry = res[which.min(res[, 2]), 1], :
## The response has five or fewer unique values. Are you sure you want to do
## regression?##
##
## GBM STEP - version 2.9
##
## Performing cross-validation optimisation of a boosted regression tree model
## for NA and using a family of bernoulli
## Using 144 observations and 6 predictors
## creating 10 initial models of 50 trees
##
## folds are stratified by prevalence
## total mean deviance = 1.3863
## tolerance is fixed at 0.0014
## ntrees resid. dev.
## 50 1.0832
## now adding trees...
## 100 0.9492
## 150 0.8849
## 200 0.8596
## 250 0.8469
## 300 0.8412
## 350 0.837
## 400 0.8405
## 450 0.843
## 500 0.8478
## 550 0.8511
## 600 0.8522
## 650 0.8579
## 700 0.859
## 750 0.8618
## 800 0.8706
## 850 0.8724
## 900 0.8772
## 950 0.8824
## 1000 0.8835
## 1050 0.8862## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Warning: glm.fit: algorithm did not converge## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred##
## mean total deviance = 1.386
## mean residual deviance = 0.694
##
## estimated cv deviance = 0.837 ; se = 0.107
##
## training data correlation = 0.771
## cv correlation = 0.714 ; se = 0.06
##
## training data AUC score = 0.937
## cv AUC score = 0.905 ; se = 0.031
##
## elapsed time - 0.02 minutes##
##
## GBM STEP - version 2.9
##
## Performing cross-validation optimisation of a boosted regression tree model
## for NA and using a family of bernoulli
## Using 144 observations and 6 predictors
## creating 10 initial models of 50 trees
##
## folds are stratified by prevalence
## total mean deviance = 1.3863
## tolerance is fixed at 0.0014
## ntrees resid. dev.
## 50 1.1361
## now adding trees...
## 100 1.0312
## 150 0.9742
## 200 0.9423
## 250 0.9203
## 300 0.9083
## 350 0.8951
## 400 0.8851
## 450 0.8817
## 500 0.8774
## 550 0.8726
## 600 0.8711
## 650 0.8671
## 700 0.8678
## 750 0.8651
## 800 0.865
## 850 0.8655
## 900 0.868
## 950 0.8681
## 1000 0.867
## 1050 0.8702
## 1100 0.8749
## 1150 0.8789
## 1200 0.879
## 1250 0.8821
## 1300 0.8851## Warning: glm.fit: algorithm did not converge
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred##
## mean total deviance = 1.386
## mean residual deviance = 0.555
##
## estimated cv deviance = 0.865 ; se = 0.112
##
## training data correlation = 0.831
## cv correlation = 0.702 ; se = 0.055
##
## training data AUC score = 0.967
## cv AUC score = 0.889 ; se = 0.032
##
## elapsed time - 0.03 minutes##
##
## GBM STEP - version 2.9
##
## Performing cross-validation optimisation of a boosted regression tree model
## for NA and using a family of bernoulli
## Using 144 observations and 6 predictors
## creating 10 initial models of 50 trees
##
## folds are stratified by prevalence
## total mean deviance = 1.3863
## tolerance is fixed at 0.0014
## ntrees resid. dev.
## 50 1.1222
## now adding trees...
## 100 0.99
## 150 0.9303
## 200 0.8975
## 250 0.8747
## 300 0.8529
## 350 0.8313
## 400 0.8167
## 450 0.8006
## 500 0.7882
## 550 0.7797
## 600 0.7724
## 650 0.7629
## 700 0.7529
## 750 0.7462
## 800 0.7396
## 850 0.7309
## 900 0.7265
## 950 0.7222
## 1000 0.7163
## 1050 0.7105
## 1100 0.7059
## 1150 0.7008
## 1200 0.7004
## 1250 0.6975
## 1300 0.6923
## 1350 0.6874
## 1400 0.6842
## 1450 0.68
## 1500 0.6766
## 1550 0.676
## 1600 0.6716
## 1650 0.6695
## 1700 0.6662
## 1750 0.6637
## 1800 0.6613
## 1850 0.6576
## 1900 0.6592
## 1950 0.6565
## 2000 0.6531
## 2050 0.6502
## 2100 0.6481
## 2150 0.6454
## 2200 0.6434
## 2250 0.6409
## 2300 0.637
## 2350 0.635
## 2400 0.6376
## 2450 0.6358
## 2500 0.6354
## 2550 0.6339
## 2600 0.6351
## 2650 0.6338
## 2700 0.6332
## 2750 0.6339
## 2800 0.6362
## 2850 0.635
## 2900 0.636
## 2950 0.6364
## 3000 0.638
## 3050 0.6391
## 3100 0.6387
## 3150 0.6392## Warning: glm.fit: algorithm did not converge
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Warning: glm.fit: algorithm did not converge## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred##
## mean total deviance = 1.386
## mean residual deviance = 0.243
##
## estimated cv deviance = 0.633 ; se = 0.09
##
## training data correlation = 0.954
## cv correlation = 0.799 ; se = 0.036
##
## training data AUC score = 0.996
## cv AUC score = 0.941 ; se = 0.017
##
## elapsed time - 0.07 minutesresults <- list.files("projection/projtest", pattern = ".tif$", recursive = T, full.names = T)
length(results)#oito algoritmos, 5 projeções (present e as demais, nem vou plotar)## [1] 160Final model projection
The final_model() function has a parameter proj_dir that allows it to be run using any projection. We can project the models fit to the present.
args(final_model)## function (species_name, algorithms = NULL, select_partitions = TRUE,
## select_par = "TSSmax", select_par_val = 0.7, weight_par = NULL,
## cut_level = c("spec_sens"), scale_models = TRUE, consensus_level = 0.5,
## models_dir = "./models", final_dir = "final_models", proj_dir = "present",
## which_models = c("raw_mean"), uncertainty = FALSE, write_final = TRUE,
## ...)
## NULLfinal_model_present <- final_model(species_name = especies[1],
models_dir = "projection/projtest",
select_par = "TSSmax",
select_par_val = 0,
proj_dir = "present",
overwrite = TRUE)## [1] "Tue Nov 12 19:07:35 2019"
## Abarema_langsdorffii
## Reading evaluation files for Abarema_langsdorffii in present
## Extracting data for Abarema_langsdorffii bioclim
## Reading models from .tif files
## selecting partitions for Abarema_langsdorffii bioclim
## Standardizing models from 0 to 1
## selected final models for Abarema_langsdorffii bioclim DONE## Extracting data for Abarema_langsdorffii brt
## Reading models from .tif files
## selecting partitions for Abarema_langsdorffii brt
## Standardizing models from 0 to 1
## selected final models for Abarema_langsdorffii brt DONE## Extracting data for Abarema_langsdorffii domain
## Reading models from .tif files
## selecting partitions for Abarema_langsdorffii domain
## Standardizing models from 0 to 1
## selected final models for Abarema_langsdorffii domain DONE## Extracting data for Abarema_langsdorffii glm
## Reading models from .tif files
## selecting partitions for Abarema_langsdorffii glm
## Standardizing models from 0 to 1
## selected final models for Abarema_langsdorffii glm DONE## Extracting data for Abarema_langsdorffii maxent
## Reading models from .tif files
## selecting partitions for Abarema_langsdorffii maxent
## Standardizing models from 0 to 1
## selected final models for Abarema_langsdorffii maxent DONE## Extracting data for Abarema_langsdorffii rf
## Reading models from .tif files
## selecting partitions for Abarema_langsdorffii rf
## Standardizing models from 0 to 1
## selected final models for Abarema_langsdorffii rf DONE## Extracting data for Abarema_langsdorffii svme
## Reading models from .tif files
## selecting partitions for Abarema_langsdorffii svme
## Standardizing models from 0 to 1
## selected final models for Abarema_langsdorffii svme DONE## Extracting data for Abarema_langsdorffii svmk
## Reading models from .tif files
## selecting partitions for Abarema_langsdorffii svmk
## Standardizing models from 0 to 1
## selected final models for Abarema_langsdorffii svmk DONE## [1] "DONE svmk !"or the projected models,
final_model_proj1 <- final_model(species_name = especies[1],
models_dir = "projection/projtest",
select_par = "TSSmax",
select_par_val = 0,
proj_dir = "proj1", overwrite = TRUE)## [1] "Tue Nov 12 19:08:10 2019"
## Abarema_langsdorffii
## Reading evaluation files for Abarema_langsdorffii in proj1
## Extracting data for Abarema_langsdorffii bioclim
## Reading models from .tif files
## selecting partitions for Abarema_langsdorffii bioclim
## Standardizing models from 0 to 1
## selected final models for Abarema_langsdorffii bioclim DONE## Extracting data for Abarema_langsdorffii brt
## Reading models from .tif files
## selecting partitions for Abarema_langsdorffii brt
## Standardizing models from 0 to 1
## selected final models for Abarema_langsdorffii brt DONE## Extracting data for Abarema_langsdorffii domain
## Reading models from .tif files
## selecting partitions for Abarema_langsdorffii domain
## Standardizing models from 0 to 1
## selected final models for Abarema_langsdorffii domain DONE## Extracting data for Abarema_langsdorffii glm
## Reading models from .tif files
## selecting partitions for Abarema_langsdorffii glm
## Standardizing models from 0 to 1
## selected final models for Abarema_langsdorffii glm DONE## Extracting data for Abarema_langsdorffii maxent
## Reading models from .tif files
## selecting partitions for Abarema_langsdorffii maxent
## Standardizing models from 0 to 1
## selected final models for Abarema_langsdorffii maxent DONE## Extracting data for Abarema_langsdorffii rf
## Reading models from .tif files
## selecting partitions for Abarema_langsdorffii rf
## Standardizing models from 0 to 1
## selected final models for Abarema_langsdorffii rf DONE## Extracting data for Abarema_langsdorffii svme
## Reading models from .tif files
## selecting partitions for Abarema_langsdorffii svme
## Standardizing models from 0 to 1
## selected final models for Abarema_langsdorffii svme DONE## Extracting data for Abarema_langsdorffii svmk
## Reading models from .tif files
## selecting partitions for Abarema_langsdorffii svmk
## Standardizing models from 0 to 1
## selected final models for Abarema_langsdorffii svmk DONE## [1] "DONE svmk !"final_model_proj2 <- final_model(species_name = especies[1],
models_dir = "projection/projtest",
select_par = "TSSmax",
select_par_val = 0,
proj_dir = "proj2", overwrite = TRUE)## [1] "Tue Nov 12 19:08:27 2019"
## Abarema_langsdorffii
## Reading evaluation files for Abarema_langsdorffii in proj2
## Extracting data for Abarema_langsdorffii bioclim
## Reading models from .tif files
## selecting partitions for Abarema_langsdorffii bioclim
## Standardizing models from 0 to 1
## selected final models for Abarema_langsdorffii bioclim DONE## Extracting data for Abarema_langsdorffii brt
## Reading models from .tif files
## selecting partitions for Abarema_langsdorffii brt
## Standardizing models from 0 to 1
## selected final models for Abarema_langsdorffii brt DONE## Extracting data for Abarema_langsdorffii domain
## Reading models from .tif files
## selecting partitions for Abarema_langsdorffii domain
## Standardizing models from 0 to 1
## selected final models for Abarema_langsdorffii domain DONE## Extracting data for Abarema_langsdorffii glm
## Reading models from .tif files
## selecting partitions for Abarema_langsdorffii glm
## Standardizing models from 0 to 1
## selected final models for Abarema_langsdorffii glm DONE## Extracting data for Abarema_langsdorffii maxent
## Reading models from .tif files
## selecting partitions for Abarema_langsdorffii maxent
## Standardizing models from 0 to 1
## selected final models for Abarema_langsdorffii maxent DONE## Extracting data for Abarema_langsdorffii rf
## Reading models from .tif files
## selecting partitions for Abarema_langsdorffii rf
## Standardizing models from 0 to 1
## selected final models for Abarema_langsdorffii rf DONE## Extracting data for Abarema_langsdorffii svme
## Reading models from .tif files
## selecting partitions for Abarema_langsdorffii svme
## Standardizing models from 0 to 1
## selected final models for Abarema_langsdorffii svme DONE## Extracting data for Abarema_langsdorffii svmk
## Reading models from .tif files
## selecting partitions for Abarema_langsdorffii svmk
## Standardizing models from 0 to 1
## selected final models for Abarema_langsdorffii svmk DONE## [1] "DONE svmk !"final_model_proj3 <- final_model(species_name = especies[1],
models_dir = "projection/projtest",
select_par = "TSSmax",
select_par_val = 0,
proj_dir = "proj3", overwrite = TRUE)## [1] "Tue Nov 12 19:09:02 2019"
## Abarema_langsdorffii
## Reading evaluation files for Abarema_langsdorffii in proj3
## Extracting data for Abarema_langsdorffii bioclim
## Reading models from .tif files
## selecting partitions for Abarema_langsdorffii bioclim
## Standardizing models from 0 to 1
## selected final models for Abarema_langsdorffii bioclim DONE## Extracting data for Abarema_langsdorffii brt
## Reading models from .tif files
## selecting partitions for Abarema_langsdorffii brt
## Standardizing models from 0 to 1
## selected final models for Abarema_langsdorffii brt DONE## Extracting data for Abarema_langsdorffii domain
## Reading models from .tif files
## selecting partitions for Abarema_langsdorffii domain
## Standardizing models from 0 to 1
## selected final models for Abarema_langsdorffii domain DONE## Extracting data for Abarema_langsdorffii glm
## Reading models from .tif files
## selecting partitions for Abarema_langsdorffii glm
## Standardizing models from 0 to 1
## selected final models for Abarema_langsdorffii glm DONE## Extracting data for Abarema_langsdorffii maxent
## Reading models from .tif files
## selecting partitions for Abarema_langsdorffii maxent
## Standardizing models from 0 to 1
## selected final models for Abarema_langsdorffii maxent DONE## Extracting data for Abarema_langsdorffii rf
## Reading models from .tif files
## selecting partitions for Abarema_langsdorffii rf
## Standardizing models from 0 to 1
## selected final models for Abarema_langsdorffii rf DONE## Extracting data for Abarema_langsdorffii svme
## Reading models from .tif files
## selecting partitions for Abarema_langsdorffii svme
## Standardizing models from 0 to 1
## selected final models for Abarema_langsdorffii svme DONE## Extracting data for Abarema_langsdorffii svmk
## Reading models from .tif files
## selecting partitions for Abarema_langsdorffii svmk
## Standardizing models from 0 to 1
## selected final models for Abarema_langsdorffii svmk DONE## [1] "DONE svmk !"final_model_proj4 <- final_model(species_name = especies[1],
models_dir = "projection/projtest",
select_par = "TSSmax",
select_par_val = 0,
proj_dir = "proj4", overwrite = TRUE)## [1] "Tue Nov 12 19:09:38 2019"
## Abarema_langsdorffii
## Reading evaluation files for Abarema_langsdorffii in proj4
## Extracting data for Abarema_langsdorffii bioclim
## Reading models from .tif files
## selecting partitions for Abarema_langsdorffii bioclim
## Standardizing models from 0 to 1
## selected final models for Abarema_langsdorffii bioclim DONE## Extracting data for Abarema_langsdorffii brt
## Reading models from .tif files
## selecting partitions for Abarema_langsdorffii brt
## Standardizing models from 0 to 1
## selected final models for Abarema_langsdorffii brt DONE## Extracting data for Abarema_langsdorffii domain
## Reading models from .tif files
## selecting partitions for Abarema_langsdorffii domain
## Standardizing models from 0 to 1
## selected final models for Abarema_langsdorffii domain DONE## Extracting data for Abarema_langsdorffii glm
## Reading models from .tif files
## selecting partitions for Abarema_langsdorffii glm
## Standardizing models from 0 to 1
## selected final models for Abarema_langsdorffii glm DONE## Extracting data for Abarema_langsdorffii maxent
## Reading models from .tif files
## selecting partitions for Abarema_langsdorffii maxent
## Standardizing models from 0 to 1
## selected final models for Abarema_langsdorffii maxent DONE## Extracting data for Abarema_langsdorffii rf
## Reading models from .tif files
## selecting partitions for Abarema_langsdorffii rf
## Standardizing models from 0 to 1
## selected final models for Abarema_langsdorffii rf DONE## Extracting data for Abarema_langsdorffii svme
## Reading models from .tif files
## selecting partitions for Abarema_langsdorffii svme
## Standardizing models from 0 to 1
## selected final models for Abarema_langsdorffii svme DONE## Extracting data for Abarema_langsdorffii svmk
## Reading models from .tif files
## selecting partitions for Abarema_langsdorffii svmk
## Standardizing models from 0 to 1
## selected final models for Abarema_langsdorffii svmk DONE## [1] "DONE svmk !"