Code owners
Assign users and groups as approvers for specific file changes. Learn more.
03_01_modelling_ssdm.R 9.13 KiB
library(dplyr)
library(tidyr)
library(sf)
library(caret)
library(cito)
library(pROC)
library(parallel)
library(doParallel)
library(foreach)
source("R/utils.R")
load("data/r_objects/model_data.RData")
species_processed = list.files("data/r_objects/ssdm_results/performance/", pattern = "RData") %>%
stringr::str_remove(".RData")
data_split = model_data %>%
sf::st_drop_geometry() %>%
dplyr::filter(!species %in% species_processed) %>%
dplyr::group_by(species) %>%
dplyr::group_split()
# ----------------------------------------------------------------------#
# Set up parallel processing ####
# ----------------------------------------------------------------------#
num_cores <- 2 # Some packages use internal parallelization, so be careful with increasing this param
cl <- makeCluster(num_cores)
registerDoParallel(cl)
# ----------------------------------------------------------------------#
# Run training ####
# ----------------------------------------------------------------------#
# Use foreach for parallel processing across species
foreach(pa_spec = data_split, .packages = c("dplyr", "caret", "cito", "pROC", "tidyr")) %dopar% {
na_performance = list(auc = NA_real_, accuracy = NA_real_, kappa = NA_real_,
precision = NA_real_, recall = NA_real_, f1 = NA_real_,
tp = NA_real_, fp = NA_real_, tn = NA_real_, fn = NA_real_)
predictors = paste0("bio", 1:19)
species = pa_spec$species[1]
print(species)
pa_spec$present_fct = factor(pa_spec$present, levels = c("0", "1"), labels = c("A", "P"))
folds = sort(unique(na.omit(pa_spec$fold_global)))
# Outer CV loop (for averaging performance metrics)
performance_cv = lapply(folds, function(fold){
print(paste("Fold", fold))
## Preparations #####
data_train = dplyr::filter(pa_spec, record_type == "background" | fold_global != fold) # include background samples
data_test = dplyr::filter(pa_spec, fold_global == fold)
if(all(data_train$present == 0)){ # If only absences in training dataset --> return
return(NULL)
}
train_ctrl = caret::trainControl(
search = "random",
classProbs = TRUE,
number = 5,
summaryFunction = caret::twoClassSummary,
savePredictions = "final"
)
## Random Forest #####
rf_performance = tryCatch({
# Fit model rf_fit = caret::train(
x = data_train[, predictors],
y = data_train$present_fct,
method = "ranger",
trControl = train_ctrl,
tuneLength = 8,
weights = data_train$weight,
verbose = F
)
evaluate_model(rf_fit, data_test)
}, error = function(e){
na_performance
})
## Gradient Boosted Machine ####
gbm_performance = tryCatch({
gbm_fit = train(
x = data_train[, predictors],
y = data_train$present_fct,
method = "gbm",
trControl = train_ctrl,
tuneLength = 8,
weights = data_train$weight,
verbose = F
)
evaluate_model(gbm_fit, data_test)
}, error = function(e){
na_performance
})
## Generalized Additive Model ####
gam_performance = tryCatch({
gam_fit = train(
x = data_train[, predictors],
y = data_train$present_fct,
method = "gamSpline",
tuneLength = 8,
weights = data_train$weight,
trControl = train_ctrl
)
evaluate_model(gam_fit, data_test)
}, error = function(e){
na_performance
})
## Neural Network ####
formula = as.formula(paste0("present ~ ", paste(predictors, collapse = '+')))
nn_performance = tryCatch({
nn_fit = dnn(
formula,
data = data_train,
hidden = c(200L, 200L, 200L),
loss = "binomial",
epochs = 200L,
burnin = 200L,
early_stopping = 25,
lr = 0.001,
batchsize = min(ceiling(nrow(data_train)/10), 64),
dropout = 0.25,
optimizer = config_optimizer("adam"),
validation = 0.2,
device = "cuda",
verbose = F,
plot = F
)
evaluate_model(nn_fit, data_test)
}, error = function(e){
na_performance })
# Summarize results
performance_summary = tibble(
species = !!species,
obs = nrow(data_train),
fold_global = fold,
model = c("rf", "gbm", "gam", "nn"),
auc = c(rf_performance$auc, gbm_performance$auc, gam_performance$auc, nn_performance$auc),
accuracy = c(rf_performance$accuracy, gbm_performance$accuracy, gam_performance$accuracy, nn_performance$accuracy),
kappa = c(rf_performance$kappa, gbm_performance$kappa, gam_performance$kappa, nn_performance$kappa),
precision = c(rf_performance$precision, gbm_performance$precision, gam_performance$precision, nn_performance$precision),
recall = c(rf_performance$recall, gbm_performance$recall, gam_performance$recall, nn_performance$recall),
f1 = c(rf_performance$f1, gbm_performance$f1, gam_performance$f1, nn_performance$f1),
tp = c(rf_performance$tp, gbm_performance$tp, gam_performance$tp, nn_performance$tp),
fp = c(rf_performance$fp, gbm_performance$fp, gam_performance$fp, nn_performance$fp),
tn = c(rf_performance$tn, gbm_performance$tn, gam_performance$tn, nn_performance$tn),
fn = c(rf_performance$fn, gbm_performance$fn, gam_performance$fn, nn_performance$fn)
) %>%
tidyr::pivot_longer(-any_of(c("species", "obs", "fold_global", "model")), names_to = "metric", values_to = "value")
return(performance_summary)
})
# Combine and save evaluation results
performance_spec = bind_rows(performance_cv)
if(nrow(performance_spec) != 0){
performance_spec = performance_spec %>%
dplyr::arrange(fold_global, model, metric)
save(performance_spec, file = paste0("data/r_objects/ssdm_results/performance/", species, ".RData"))
}
}
# Stop the parallel cluster when done
stopCluster(cl)
# Combine results
files = list.files("data/r_objects/ssdm_results/performance/", full.names = T, pattern = ".RData")
ssdm_performance = lapply(files, function(f){load(f); return(performance_spec)}) %>%
bind_rows()
save(ssdm_performance, file = "data/r_objects/ssdm_performance.RData")
# ----------------------------------------------------------------------#
# Train full models ####
# ----------------------------------------------------------------------#
data_split = model_data %>%
dplyr::mutate(present_fct = factor(present, levels = c("0", "1"), labels = c("A", "P"))) %>%
dplyr::group_by(species) %>%
dplyr::group_split()
for(pa_spec in data_split){
## Preparations ####
species = pa_spec$species[1]
print(as.character(species))
# Create inner CV folds for model training
cv_train = blockCV::cv_spatial(
pa_spec,
column = "present",
k = 5,
progress = F, plot = F, report = F
)
pa_spec$fold_train = cv_train$folds_ids
# Drop geometry
pa_spec$geometry = NULL
# Define caret training routine index_train = lapply(unique(sort(pa_spec$fold_train)), function(x){
return(which(pa_spec$fold_train != x))
})
train_ctrl = caret::trainControl(
search = "random",
classProbs = TRUE,
index = index_train,
summaryFunction = caret::twoClassSummary,
savePredictions = "final"
)
# Define predictors
predictors = paste0("bio", 1:19)
# Random Forest ####
try({
rf_fit = caret::train(
x = pa_spec[, predictors],
y = pa_spec$present_fct,
method = "ranger",
trControl = train_ctrl,
tuneLength = 8,
weights = pa_spec$weight,
verbose = F
)
save(rf_fit, file = paste0("data/r_objects/ssdm_results/full_models/", species, "_rf_fit.RData"))
})
# Gradient Boosted Machine ####
try({
gbm_fit = train(
x = pa_spec[, predictors],
y = pa_spec$present_fct,
method = "gbm",
trControl = train_ctrl,
tuneLength = 8,
weights = pa_spec$weight,
verbose = F
)
save(gbm_fit, file = paste0("data/r_objects/ssdm_results/full_models/", species, "_gbm_fit.RData"))
})
# Generalized Additive Model ####
try({
gam_fit = train(
x = pa_spec[, predictors],
y = pa_spec$present_fct,
method = "gamSpline",
tuneLength = 8,
weights = pa_spec$weight,
trControl = train_ctrl
)
save(gam_fit, file = paste0("data/r_objects/ssdm_results/full_models/", species, "_gam_fit.RData"))
})
# Neural Network ####
try({
formula = as.formula(paste0("present ~ ", paste(predictors, collapse = '+')))
nn_fit = dnn(
formula,
data = pa_spec,
hidden = c(200L, 200L, 200L),
loss = "binomial",
epochs = 200L,
burnin = 200L,
early_stopping = 25,
lr = 0.001,
batchsize = min(ceiling(nrow(pa_spec)/10), 64), dropout = 0.25,
optimizer = config_optimizer("adam"),
validation = 0.2,
device = "cuda",
verbose = F,
plot = F
)
save(nn_fit, file = paste0("data/r_objects/ssdm_results/full_models/", species, "_nn_fit.RData"))
})
}