# Get a unique list of species names to iterate through and process each one.
sp_names <- unique(psa$species)
# Loop through each species name in the 'sp_names' vector.
for (i in 1:length(sp_names)) {
# Get the current species name for this iteration.
s_name <- sp_names[i]
# Print a message to the console to track progress.
message(paste0("\n===== Processing species: ", s_name, " ====="))
## %######################################################%##
#### 1. Prepare Data for the Current Species ####
## %######################################################%##
# Load the environmental raster data specific to the species.
env_r <- terra::rast(sp_pc[s_name])
# Get the names of the environmental variables (predictors).
var_names <- names(env_r)
# Filter the presence-absence data for the current species.
sp_pa <- psa %>% dplyr::filter(species == s_name)
# Filter the background data for the current species.
sp_bg <- bkgrnd %>% dplyr::filter(species == s_name)
# Extract environmental data for presence-absence points, removing any points with NA values.
sp_pa <- sdm_extract(
sp_pa,
x = "x",
y = "y",
env_layer = env_r,
filter_na = TRUE
)
# Extract environmental data for background points.
sp_bg <- sdm_extract(
sp_bg,
x = "x",
y = "y",
env_layer = env_r,
filter_na = TRUE
)
## %######################################################%##
#### 2. Fit Models ####
## %######################################################%##
message(" > Fitting models...")
# ---- Random Forest (tune_raf) ----
# Tune hyperparameters for a Random Forest model.
m_raf <- tune_raf(
data = sp_pa, # Presence-absence data
response = "pr_ab", # The column with presence/absence (1/0)
predictors = var_names, # Names of environmental variables
partition = ".part", # Column for data partitioning (training/testing)
grid = expand.grid(
# Define the grid of hyperparameters to test
mtry = seq(1, length(var_names), 1), # Number of variables to sample at each split
ntree = c(400, 600, 800, 1000) # Number of trees to grow
),
thr = "max_sorensen", # Threshold to convert continuous prediction to binary
metric = "SORENSEN", # Metric to evaluate the model performance
n_cores = 4 # Number of CPU cores for parallel processing
)
# ---- Maxent (tune_max) ----
# Tune hyperparameters for a Maxent model.
m_max <- tune_max(
data = sp_pa,
response = "pr_ab",
predictors = var_names,
background = sp_bg, # Background data is required for Maxent
partition = ".part",
grid = expand.grid(
# Define the grid of hyperparameters to test
regmult = seq(0.1, 2, 0.5), # Regularization multiplier to prevent overfitting
classes = c("lq", "lqh", "lqhp", "lqhpt") # Feature classes (l=linear, q=quadratic, etc.)
),
thr = "max_sorensen",
metric = "SORENSEN",
n_cores = 4
)
# ---- Boosted Regression Tree (tune_gbm) ----
# Tune hyperparameters for a Boosted Regression Tree model.
m_gbm <- tune_gbm(
data = sp_pa,
response = "pr_ab",
predictors = var_names,
partition = ".part",
grid = expand.grid(
# Define the grid of hyperparameters to test
n.trees = seq(150, 200, 10), # Number of trees
shrinkage = seq(1, 1.5, 0.2), # Learning rate
n.minobsinnode = seq(1, 15, 2) # Minimum number of observations in a terminal node
),
thr = "max_sorensen",
metric = "SORENSEN",
n_cores = 4
)
# ---- Generalized Additive Model (fit_gam) ----
# Determine the number of training samples.
n_t <- flexsdm:::n_training(data = sp_pa, partition = ".part")
# Set a starting number for basis functions (k).
candidate_k <- 20
# Reduce 'k' until it's smaller than the number of training samples to avoid errors.
while (
any(
n_t <
flexsdm:::n_coefficients(
data = sp_pa,
predictors = var_names,
k = candidate_k
)
)
) {
candidate_k <- candidate_k - 3
}
# Fit a GAM model with the adjusted 'k'.
m_gam <- fit_gam(
data = sp_pa,
response = "pr_ab",
predictors = var_names,
partition = ".part",
thr = "max_sorensen",
k = candidate_k # Basis functions parameter
)
# ---- Generalized Linear Model (fit_glm) ----
# Fit a GLM with a second-degree polynomial term.
m_glm <- fit_glm(
data = sp_pa,
response = "pr_ab",
predictors = var_names,
partition = ".part",
thr = "max_sorensen",
poly = 2 # Degree of polynomial functions
)
# ---- Gaussian Process (fit_gau) ----
# Fit a Gaussian Process model.
m_gau <- fit_gau(
data = sp_pa,
response = "pr_ab",
predictors = var_names,
partition = ".part",
thr = "max_sorensen"
)
# ---- Neural Network (tune_net) ----
# Tune hyperparameters for a Neural Network model.
m_net <- tune_net(
data = sp_pa,
response = "pr_ab",
predictors = var_names,
partition = ".part",
grid = expand.grid(
# Define the grid of hyperparameters to test
size = 2:length(var_names), # Number of units in the hidden layer
decay = c(seq(0.01, 1, 0.05), 1, 3, 4, 5, 6) # Weight decay for regularization
),
thr = "max_sorensen",
metric = "SORENSEN",
n_cores = 4
)
# ---- Support Vector Machine (tune_svm) ----
# Tune hyperparameters for a Support Vector Machine model.
m_svm <- tune_svm(
data = sp_pa,
response = "pr_ab",
predictors = var_names,
partition = ".part",
grid = expand.grid(
C = seq(2, 60, 5), # Cost parameter for regularization
sigma = c(seq(0.001, 0.3, 0.01)) # Kernel parameter
),
thr = "max_sorensen",
metric = "SORENSEN",
n_cores = 4
)
# Gather all fitted model objects (those starting with "m_") into a single named list.
model_obj_list <- mget(grep("^m_", ls(), value = TRUE))
# Clean up the names in the list (e.g., "m_raf" becomes "raf").
names(model_obj_list) <- gsub("m_", "", names(model_obj_list))
## %######################################################%##
#### 3. Partial Dependence & Variable Importance ####
## %######################################################%##
message(" > Calculating PDP and Variable Importance...")
# Initialize an empty list to store variable importance results for each model.
varimp <- list()
# Loop through each fitted model to generate plots and calculate importance.
for (i in 1:length(model_obj_list)) {
mod_name <- names(model_obj_list)[i]
model_obj <- model_obj_list[[i]]
# Generate Partial Dependence Plots (PDP) to visualize variable effects.
# The 'p_pdp' function extracts the core model (e.g., randomForest) from the flexsdm object.
p <- p_pdp(model = model_obj$model, training_data = sp_pa) +
ggtitle(paste(s_name, "-", mod_name)) # Add a title to the plot
# If the plot was created successfully, save it as a PNG file.
if (!is.null(p)) {
ggsave(
plot = p,
filename = file.path(
dir_fig,
paste0(s_name, "_", mod_name, "_pdp.png")
),
width = 8,
height = 6,
dpi = 300
)
}
# Calculate variable importance using a permutation-based method.
varimp[[i]] <- flexsdm::sdm_varimp(
model = model_obj,
data = sp_pa,
response = "pr_ab",
predictors = var_names,
n_sim = 30, # Number of permutations to run
n_cores = 2,
thr = "max_sorensen",
clamp = TRUE,
pred_type = "cloglog" # Use cloglog output for predictions
)
}
# Combine the variable importance results from all models into a single data frame.
varimp <- dplyr::bind_rows(varimp)
# Save the combined variable importance table to a text file.
readr::write_tsv(varimp, file.path(dir_perf, paste0(s_name, "_varimp.txt")))
# Create a bar plot to visualize variable importance across models.
varimp2 <- varimp %>%
tidyr::pivot_longer(
cols = TPR:IMAE,
names_to = "metric",
values_to = "value"
) %>%
filter(metric %in% c("AUC", "TSS", "SORENSEN")) # Filter for specific metrics
p_vi <- ggplot(varimp2, aes(x = predictors, y = value, fill = model)) +
scale_fill_brewer(palette = "Set1") +
geom_bar(stat = "identity", position = position_dodge()) +
facet_grid(model ~ metric, scales = "free") +
coord_flip() +
labs(
x = "Predictor variables",
y = "Variable importance"
) +
theme_bw()
print(p_vi)
ggsave(
filename = file.path(dir_fig, paste0(s_name, "_varImp.png")),
plot = p_vi,
width = 12,
height = 14,
dpi = 300,
unit = "cm",
scale = 1.5
)
# Summarize performance metrics (e.g., SORENSEN, TSS) for all individual models.
model_perf <- flexsdm::sdm_summarize(model_obj_list)
# Save the performance summary to a text file.
readr::write_tsv(
model_perf,
file.path(dir_perf, paste0(s_name, "_models_performance.txt"))
)
## %######################################################%##
#### 4. Filter Models for Ensemble ####
## %######################################################%##
message(" > Filtering models to perform ensemble...")
# Filter models based on performance criteria to select the "best" ones for the ensemble.
# Here, we select models with a mean Sorensen score >= 0.7.
# We also exclude models with thresholds of 0 or 1, as they might indicate overfitting.
good_models_names <- model_perf$model[
model_perf$SORENSEN_mean >= 0.7 &
model_perf$thr_value > 0 &
model_perf$thr_value < 1
]
# Create a new list containing only the selected "good" models.
good_models <- model_obj_list[good_models_names]
## %######################################################%##
#### 5. Ensemble Modeling ####
## %######################################################%##
# Initialize an empty list to store ensemble models.
ensemble_list <- list()
# Only proceed if there is more than one "good" model to create an ensemble.
if (length(good_models) > 1) {
message(" > Fitting model ensemble...")
# Create an ensemble by averaging the predictions of the good models.
m_mean <- flexsdm::fit_ensemble(
good_models,
ens_method = "mean", # Ensemble method: mean, weighted_mean, median, etc.
thr_model = "max_sorensen", # Thresholding method for individual models
thr = "max_sorensen", # Thresholding method for the final ensemble
metric = "SORENSEN"
)
# Create an ensemble using the median of the predictions.
m_median <-
flexsdm::fit_ensemble(
good_models,
ens_method = "median",
thr_model = "max_sorensen",
thr = "max_sorensen",
metric = "SORENSEN"
)
# Combine the ensemble models into a list.
ensemble_list <- list(mean = m_mean, median = m_median)
# Summarize the performance of the ensemble models.
ens_perf <- flexsdm::sdm_summarize(ensemble_list)
# Combine the performance tables of individual models and ensemble models.
model_perf_with_ens <- dplyr::bind_rows(model_perf, ens_perf)
# Overwrite the previous performance file to include the ensemble results.
readr::write_tsv(
model_perf_with_ens,
file.path(dir_perf, paste0(s_name, "_models_performance.txt"))
)
}
## %######################################################%##
#### 6. Spatial Predictions ####
## %######################################################%##
message(" > Predicting models ...")
names(good_models) <- NULL
# Predict habitat suitability for each of the "good" individual models.
if (length(good_models) > 0) {
prd_list <- flexsdm::sdm_predict(
models = good_models, # List of good models
pred = env_r, # Environmental raster layers to predict onto
thr = "max_sorensen",
clamp = TRUE, # Restrict predictions to the range of training data
con_thr = TRUE, # Produce continuous and binary (thresholded) predictions
pred_type = "cloglog"
)
# Loop through the list of predictions and save each one as a GeoTIFF raster file.
for (alg in names(prd_list)) {
terra::writeRaster(
prd_list[[alg]],
filename = file.path(dir_current, alg, paste0(s_name, ".tif")),
overwrite = TRUE,
filetype = "GTiff"
)
}
}
# Predict habitat suitability for the ensemble models.
if (length(ensemble_list) > 0) {
for (i in seq_along(ensemble_list)) {
ens_prd <- sdm_predict(
models = ensemble_list[[i]],
pred = env_r,
thr = "max_sorensen",
con_thr = TRUE,
clamp = TRUE,
pred_type = "cloglog"
)
# Save prediction as a GeoTIFF.
writeRaster(
ens_prd[[1]],
filename = file.path(
dir_current,
names(ensemble_list[i]),
paste0(s_name, ".tif")
),
overwrite = TRUE,
filetype = "GTiff"
)
}
}
## %######################################################%##
#### 7. Future projections ####
## %######################################################%##
message(" > Predicting models for future projections ...")
fut_env <- fut_sp_pc_list[[s_name]]
# Loop over each scenario raster stack (fut_sp_pc_list already has species->scenario->raster)
for (sc_name in names(fut_env)) {
env_r <- fut_env[[sc_name]] # raster stack for this scenario
# Create output folder if it doesn't exist
dir_future <- proj_dir[grep(paste0("2_Outputs/2_Projection"), proj_dir)][1]
# --- Predict individual models ---
if (length(good_models) > 0) {
prd_list <- flexsdm::sdm_predict(
models = good_models,
pred = env_r,
thr = "max_sorensen",
clamp = TRUE,
con_thr = TRUE,
pred_type = "cloglog"
)
# Save individual model predictions
for (alg in names(prd_list)) {
terra::writeRaster(
prd_list[[alg]],
file.path(
dir_future,
sc_name,
paste0("Algorithm"),
alg,
paste0(s_name, ".tif")
),
overwrite = TRUE,
filetype = "GTiff"
)
}
}
# --- Predict ensemble models ---
if (length(ensemble_list) > 0) {
for (i in seq_along(ensemble_list)) {
ens_prd <- flexsdm::sdm_predict(
models = ensemble_list[[i]],
pred = env_r,
thr = "max_sorensen",
clamp = TRUE,
con_thr = TRUE,
pred_type = "cloglog"
)
# Save ensemble prediction
writeRaster(
ens_prd[[1]],
filename = file.path(
dir_future,
sc_name,
paste0("Algorithm"),
paste0(names(ensemble_list[i])),
paste0(s_name, ".tif")
),
overwrite = TRUE,
filetype = "GTiff"
)
}
}
}
## %######################################################%##
#### 8. Clean Up for Next Iteration ####
## %######################################################%##
# Remove all model objects (starting with "m_") to free up memory.
rm(list = grep("^m_", ls(), value = TRUE))
# Remove another temporary object if it exists.
if (exists("alg_prd_cont")) {
rm(alg_prd_cont)
}
# Trigger garbage collection to release memory explicitly.
gc()
} # End of the main species loop.