# 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(embarcadero)  # building and analysing BART models
library(fuzzySim)  # computing favourability

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


# IMPORT DATA ####

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

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

layers <- terra::rast(paste0("../outputs/", my_species, "/layers_cut/layers_cut.tif"))
plot(layers[[1:4]])

dat <- read.csv(paste0("../outputs/", my_species, "/occ_data/presences_gridded.csv"))
nrow(dat)
head(dat)
names(dat)

# convert to SpatVector and map:
dat_sv <- terra::vect(dat, geom = c("x", "y"), keepgeom = TRUE, crs = "EPSG:4326")
plot(dat_sv, "presence", cex = 0.5)
plot(subset(dat_sv, dat_sv$selected == TRUE), "presence", cex = 0.5)

# define the modelling columns:
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

table(dat[ , spc_col])  # number of presences and absences
table(dat[dat$selected, spc_col])  # number of presences and selected absences

# check for (and remove) variables with no variation in the study region:

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

constants

if (length(constants) > 0) var_cols <- var_cols[-constants]


# COMPUTE FULL MODEL (all variables) ####

# we'll compute models with selected variables later but, for now, let's build one with all our candidate variables:

mod_full <- dbarts::bart(x.train = dat[dat$selected, var_cols], y.train = dat[dat$selected, spc_col], keeptrees = TRUE, seed = 123)

summary(mod_full)


# save the model object to disk:

invisible(mod_full$fit$state)  # if you want to use this BART model e.g. for prediction or for plotting response curves IN FUTURE R SESSIONS, you need to run this command to explicitly ask for the full information to be included when you save the model object - see "Saving" section under ?bart Details

dir.create(paste0("../outputs/", my_species, "/models/"), showWarnings = FALSE)
saveRDS(mod_full, paste0("../outputs/", my_species, "/models/mod_full.rds"))


# OPTIONALLY, FINE-TUNE THE PRIOR PARAMETERS ####

# BART models are good at providing adequate priors with the default settings, but you may achieve slightly higher predictive performance if you 'retune' the prior parameters via cross-validation (mind this may take a VERY long time to run!):

# set.seed(1234)
# mod_full_retuned <- retune(mod_full)
# summary(mod_full_retuned)

# invisible(mod_full_retuned$fit$state)
# saveRDS(mod_full_retuned, paste0("../outputs/", my_species, "/models/mod_full_retuned.rds"))


# GET PREDICTIONS ON THE DATA TABLE ####

source("https://raw.githubusercontent.com/AMBarbosa/unpackaged/master/predict_bart_df")  # I edited the 'embarcadero' predict function to work on data frames and not only raster layers

preds <- predict_bart_df(mod_full, dat, quantiles = c(0.05, 0.95))
# takes time for large datasets!

head(preds)  # mean posterior prediction, and the lower and upper bounds of the credibility interval specified in 'quantiles' above

preds$uncert <- preds[ , 3] - preds[ , 2]  # width of the credibility interval = uncertainty of the prediction at each site

# change to more self-explanatory column names:
names(preds) <- c("BART_P", "BART_P_lower", "BART_P_upper", "BART_P_uncert")

# append BART predictions to the data frame:
dat <- data.frame(dat, preds)
head(dat)


# convert to favourability (presence probability without the effect of modelled prevalence):

dat$BART_F <- fuzzySim::Fav(pred = dat$BART_P, sample.preval = fuzzySim::prevalence(model = mod_full))

dat$BART_F_lower <- fuzzySim::Fav(pred = dat$BART_P_lower, sample.preval = fuzzySim::prevalence(model = mod_full))

dat$BART_F_upper <- fuzzySim::Fav(pred = dat$BART_P_upper, sample.preval = fuzzySim::prevalence(model = mod_full))

dat$BART_F_uncert <- dat$BART_F_upper - dat$BART_F_lower  # width of the credibility interval = uncertainty of the prediction at each site

head(dat)


# save predictions to a CSV file:

head(dat)
names(dat)
pred_columns <- grep("cells|_P|_F", names(dat))
names(dat)[pred_columns]

dir.create(paste0("../outputs/", my_species, "/predictions/"), showWarnings = FALSE)
write.csv(dat[ , pred_columns], file = paste0("../outputs/", my_species, "/predictions/BART_preds.csv"), row.names = FALSE)


# map the predictions:

# convert to SpatVector (spatial object):
dat_sv <- terra::vect(dat, geom = c("x", "y"), keepgeom = TRUE, crs = "EPSG:4326")

plot(dat_sv, "BART_P", cex = 0.5, type = "continuous", range = c(0, 1), main = "Presence probability")

plot(dat_sv, "BART_F", cex = 0.5, type = "continuous", range = c(0, 1), main = "Environmental favourability")


# map the predictions with credible intervals:

par(oma = c(0, 0, 1.5, 0))
plot(dat_sv, c("BART_F", "BART_F_uncert", "BART_F_lower", "BART_F_upper"), cex = 0.5, type = "continuous", range = c(0, 1), mar = c(1.5, 1, 1.5, 1))
title(my_species, outer = TRUE)


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


# GET PREDICTIONS ON RASTER MAPS ####
# (if you have also raster layers for your variables)

layers  # "SpatRaster" class, while ?embarcadero::predict2.bart still requires the older "RasterStack" class

BART_P_rast <- embarcadero::predict2.bart(mod_full, raster::stack(layers), quantiles = c(0.05, 0.95))  # takes time for large or high-resolution rasters!

BART_P_rast <- terra::rast(BART_P_rast)  # back to modern "SpatRaster" class

names(BART_P_rast)
names(BART_P_rast) <- c("BART_P", "BART_P_lower", "BART_P_upper")

BART_P_rast$uncertainty <- BART_P_rast[[3]] - BART_P_rast[[2]]

plot(BART_P_rast, range = c(0, 1), mar = c(1.5, 1, 1.5, 1))


# convert to favourability:

BART_F_rast <- fuzzySim::Fav(pred = BART_P_rast[[1:3]], sample.preval = fuzzySim::prevalence(model = mod_full))

names(BART_F_rast)

names(BART_F_rast) <- sub("_P", "_F", names(BART_F_rast))

BART_F_rast$uncertainty <- BART_F_rast[[3]] - BART_F_rast[[2]]

par(oma = c(0, 0, 1, 0))
plot(BART_F_rast, range = c(0, 1))
title(my_species, outer = TRUE)


# save the predictions to disk:

writeRaster(BART_P_rast, paste0("../outputs/", my_species, "/predictions/BART_full_pred.tif"), gdal = c("COMPRESS=DEFLATE"), overwrite = TRUE)
writeRaster(BART_F_rast, paste0("../outputs/", my_species, "/predictions/BART_full_fav.tif"), gdal = c("COMPRESS=DEFLATE"), overwrite = TRUE)
# note you can visualize these GeoTIFs with GIS software, not with ordinary image software