# R scripts for short course "Species distribution models with Bayesian statistics in R"
# A. Marcia Barbosa (https://modtools.wordpress.com)


# SETUP ####

library(terra)  # mapping
library(modEvA)  # model evaluation
library(blockCV)  # spatial block cross-validation

setwd(dirname(rstudioapi::getActiveDocumentContext()$path))


# import data:

# (you can skip this section if you just ran the entire script 2a for the same species and region)

my_species <- "Loxodonta africana"  # or another species for which you ran the previous scripts

# example marine species:
# my_species <- "Phocoena sinus"

dat <- read.csv(paste0("../outputs/", my_species, "/occ_data/presences_gridded.csv"))
preds <- read.csv(paste0("../outputs/", my_species, "/predictions/BART_preds.csv"))
dat <- merge(dat, preds, by = "cells")

dat_sv <- terra::vect(dat, geom = c("x", "y"), keepgeom = TRUE, crs = "EPSG:4326")
plot(dat_sv, "presence", cex = 0.5)
plot(dat_sv, "BART_P", cex = 0.5, type = "continuous")

spc_col <- "presence"  # name of the species presence/absence column IN THE EXAMPLE DATASET (change as appropriate)
names(dat)
var_cols <- 5:23  # variables are in these column positions IN THE EXAMPLE DATASET - change as appropriate, ask me for help if needed!
# var_cols <- c(5:9)  # in the marine example
names(dat)[var_cols]  # check that it's OK!

constants <- which(sapply(dat[ , var_cols], function(x) length(unique(x)) <= 1))
if (length(constants) > 0) var_cols <- var_cols[-constants]

vars_rast <- terra::rast(paste0("../outputs/", my_species, "/layers_cut/layers_cut.tif"))

source("https://raw.githubusercontent.com/AMBarbosa/unpackaged/master/predict_bart_df")


# EVALUATE MODEL PREDICTIONS ####

par(mfrow = c(1, 1), mar = c(5, 4, 2, 1))


## area under the ROC Curve (discrimination performance) ####

modEvA::AUC(obs = dat[ , spc_col], pred = dat[ , "BART_P"])


## threshold-based classification metrics ####

modEvA::threshMeasures(obs = dat[ , spc_col], pred = dat[ , "BART_P"], thresh = "maxTSS", measures = c("CCR", "Sensitivity", "Specificity", "Precision", "Recall", "TSS", "kappa"))


## Miller calibration line ####

modEvA::MillerCalib(obs = dat[ , spc_col], pred = dat[ , "BART_P"])


# BUT THIS JUST EVALUATES HOW THE MODELS FIT THE SAME DATA ON WHICH THEY WERE TRAINED
# YOU CAN SET ASIDE SOME DATA TO LEAVE OUT OF MODEL TRAINING AND USE FOR TESTING OUT-OF-SAMPLE PREDICTIVE CAPACITY
# BLOCK CROSS-VALIDATION (below) IS CURRENTLY THE MOST ACCEPTED METHOD


# DIVIDE STUDY AREA INTO SPATIAL BLOCKS ####
# for model cross-validation

names(dat_sv)

# get spatial blocks of a given size:
mod_region <- terra::as.polygons(vars_rast * 0)
plot(mod_region)
block_size <- sqrt(terra::expanse(mod_region) / 20)  # for example!
blocks <- blockCV::cv_spatial(dat_sv, column = spc_col, r = vars_rast[[1]], k = 5, size = block_size, seed = 321)

blocks$folds_ids

dat$fold <- dat_sv$fold <- blocks$folds_ids
head(dat)
head(dat_sv)

# check how many data points per fold:
table(dat$fold)
# and how many presences per fold:
table(subset(dat, presence == 1, select = "fold"))

# map each fold:
plot(dat_sv, "fold", cex = 0.4, alpha = 0.5, mar = c(1.5, 1, 1.5, 2), main = my_species, col = map.pal("wave"))

# add the presences:
points(dat_sv[dat_sv$presence == 1, ], cex = 0.3, pch = 4)

add_legend("topleft", legend = table(subset(dat, presence == 1, select = "fold")), title = "presences\nper fold:", cex = 0.7)


#### ASSIGNMENT 5: if you're using different data than the example, post the map you obtained with the above commands


# COMPUTE MODELS AND GET PREDICTIONS FOR CROSS-VALIDATION ####

names(dat)

folds <- sort(unique(dat_sv$fold))

for (f in folds) {
  message("modelling outside fold ", f)  # inform of progress
  dat_train <- subset(dat, dat$fold != f)
  mod_BART_fold <- dbarts::bart(x.train = dat_train[ , var_cols], y.train = dat_train[ , spc_col], keeptrees = TRUE, verbose = FALSE, seed = 123)
  dat[ , paste0("BART_fold", f, "_P")] <- predict_bart_df(mod_BART_fold, dat)
}

# see the new predictions added to the data frame:
head(dat)


# save the folds and fold predictions to a CSV file:
head(dat)
names(dat)
preds_folds <- dat[ , grep("cells|fold", names(dat))]
head(preds_folds)
write.csv(preds_folds, file = paste0("../outputs/", my_species, "/predictions/BART_folds.csv"), row.names = FALSE)


# EVALUATE EACH MODEL ON ITS VALIDATION FOLD ####

# identify columns containing the predictions of different folds:
fold_cols <- grep("_fold", names(dat))
fold_cols
names(dat)[fold_cols]  # check that it's OK!

# choose some evaluation measures to compute:
measures <- c("AUC", "TSS", "MCS")

# create an empty table to receive the cross-validation results:
crossval <- as.data.frame(matrix(nrow = length(folds), ncol = length(measures)))
colnames(crossval) <- measures
crossval  # for now it's only filled with NAs

par(mfrow = c(3, 1), mar = c(4, 3, 1, 1), oma = c(0, 0, 2, 0))

for (f in folds) {
  fold_col <- names(dat)[grep(paste0("BART_fold", f), names(dat))]
  fold_dat <- subset(dat, dat$fold == f)
  crossval[f, "AUC"] <- modEvA::AUC(obs = fold_dat[ , spc_col], pred = fold_dat[ , fold_col], simplif = TRUE, plot = TRUE, main = "AUC")
  crossval[f, "TSS"] <- modEvA::threshMeasures(obs = fold_dat[ , spc_col], pred = fold_dat[ , fold_col], thresh = "maxTSS", measures = "TSS", simplif = TRUE, standardize = FALSE, main = "TSS", ylim = c(0, 1))
  crossval[f, "MCS"] <- modEvA::MillerCalib(obs = fold_dat[ , spc_col], pred = fold_dat[ , fold_col], main = "Miller line")$slope
  mtext(paste("Fold", f), outer = TRUE, cex = 1.5)
}
# press the back arrow on the top left of your plotting pane to see the different fold evaluations

# see the previously created table, now filled with the cross-validation results:
crossval

# boxplots of cross-validation (CV) metrics:
par(mfrow = c(1, 1), mar = c(7, 3, 2, 1))
boxplot(crossval, las = 2, main = my_species)
abline(h = 1, col = "blue", lty = 2)  # remember Miller calibration slope (MCS) should ideally be close to 1 (not bigger = better)


#### ASSIGNMENT 6: if you're using different data than the example, post the figure you obtained with the above commands