# =========================================
# SPECIES DISTRIBUTION MODELLING AND
# CLIMATIC NICHE DYNAMICS OF
# Zaprionus indianus
# =========================================
# Authors: Avirup Chakraborty & Subhash Rajpurohit
# Date: April, 2026
# =========================================

# =========================================
# SECTION 1: SETUP
# =========================================

packages <- c("rgbif", "dplyr", "CoordinateCleaner", "terra",
              "corrplot", "caret", "usdm", "maxnet",
              "pROC", "PresenceAbsence", "geodata",
              "ade4", "ecospat", "purrr", "spThin",
              "ggplot2", "gridExtra", "rasterVis", "RColorBrewer")

install.packages(setdiff(packages, rownames(installed.packages())))

library(rgbif)
library(dplyr)
library(CoordinateCleaner)
library(terra)
library(corrplot)
library(caret)
library(usdm)
library(maxnet)
library(pROC)
library(PresenceAbsence)
library(geodata)
library(ade4)
library(ecospat)
library(purrr)
library(spThin)
library(ggplot2)
library(gridExtra)

set.seed(42)

# Create output directory
dir.create("C:/Users/Avirup/OneDrive/Desktop/Zap_SDM/outputs",
           showWarnings = FALSE, recursive = TRUE)

# =========================================
# SECTION 2: OCCURRENCE DATA DOWNLOAD
# =========================================
# Using occ_download() for bulk data retrieval.
# This generates a citable DOI for the occurrence dataset,
# which is required by high-impact journals.
# Credentials are set via options() for rgbif >= 3.8

options(
  gbif_user  = "avirup27",
  gbif_pwd   = "Avirup@270100",
  gbif_email = "avirup.c@ahduni.edu.in"
)

# Submit download request
download_key <- occ_download(
  pred("taxonKey", 1557446),        # GBIF taxon key for Zaprionus indianus
  pred("hasCoordinate", TRUE),
  pred("hasGeospatialIssue", FALSE),
  format = "SIMPLE_CSV"
)

# Wait for GBIF to prepare the file
occ_download_wait(download_key)

# Download and import
occ_download_get(download_key,
                 path = "C:/Users/Avirup/OneDrive/Desktop/Zap_SDM/")

df_raw <- occ_download_import(
  occ_download_get(download_key,
                   path = "C:/Users/Avirup/OneDrive/Desktop/Zap_SDM/"))

# Save the download DOI for citation in manuscript
cat("GBIF download DOI — cite this in your Methods section:\n")
print(occ_download_meta(download_key)$doi)

# Filter and clean columns
df <- df_raw %>%
  dplyr::select(decimalLongitude, decimalLatitude, basisOfRecord) %>%
  dplyr::rename(lon = decimalLongitude, lat = decimalLatitude) %>%
  dplyr::filter(!is.na(lon), !is.na(lat), lon != 0, lat != 0) %>%
  dplyr::filter(basisOfRecord %in% c("HUMAN_OBSERVATION", "PRESERVED_SPECIMEN")) %>%
  dplyr::distinct(lon, lat, .keep_all = TRUE)

cat("Records after basic filtering:", nrow(df), "\n")

# =========================================
# SECTION 3: COORDINATE CLEANING
# =========================================
# CoordinateCleaner removes records flagged as:
# centroids of countries/capitals, GBIF headquarters,
# biodiversity institutions, sea coordinates, and zero coordinates.

df$species <- "Zaprionus_indianus"

clean <- clean_coordinates(df, lon = "lon", lat = "lat")

occ <- clean[clean$.summary == TRUE, c("lon", "lat")]
occ <- as.data.frame(occ)

cat("Occurrences after coordinate cleaning:", nrow(occ), "\n")

# =========================================
# SECTION 4: SPATIAL THINNING
# =========================================
# Spatial thinning is performed before environmental data loading
# to ensure that subsequent variable selection (correlation and VIF)
# is computed on a spatially unbiased set of occurrence points,
# free from the influence of clustered sampling effort.
# A geographic distance threshold of 10 km is used (spThin),
# which is resolution-independent and fully decoupled from
# any subsequent modelling grain decisions.

occ$species <- "Zaprionus_indianus"

thinned <- thin(
  loc.data    = occ,
  lat.col     = "lat",
  long.col    = "lon",
  spec.col    = "species",
  thin.par    = 10,
  reps        = 10,
  out.dir     = tempdir(),
  write.files = TRUE,
  verbose     = FALSE
)

# Read thinned files directly (spThin v0.2.0 compatibility)
thin_files <- list.files(
  path       = tempdir(),
  pattern    = "thinned_data_thin.*\\.csv",
  full.names = TRUE
)

thin_list <- lapply(thin_files, read.csv)
best_rep  <- thin_list[[which.max(sapply(thin_list, nrow))]]
colnames(best_rep) <- tolower(colnames(best_rep))

occ_thin <- as.data.frame(best_rep[, c("lon", "lat")])

# Clean up temp files
file.remove(thin_files)

cat("Occurrences before thinning:", nrow(occ), "\n")
cat("Occurrences after thinning:",  nrow(occ_thin), "\n")

# =========================================
# SECTION 5: ENVIRONMENTAL DATA
# =========================================

files <- list.files("C:/Users/Avirup/OneDrive/Desktop/Zap_SDM/climate/wc2.1_2.5m",
                    pattern = ".tif$", full.names = TRUE)

env <- rast(files)
names(env) <- paste0("bio", 1:nlyr(env))

cat("Environmental raster loaded:", nlyr(env), "layers\n")
cat("Resolution:", res(env), "\n")
cat("Extent:", as.character(ext(env)), "\n")

# =========================================
# SECTION 6: VARIABLE SELECTION
# =========================================
# Correlation and VIF analyses are performed on spatially thinned
# occurrence points to ensure multicollinearity assessment reflects
# an unbiased sample of the environmental space.

# Extract values at thinned occurrences
env_vals <- extract(env, as.matrix(occ_thin[, c("lon", "lat")]))[, names(env)]
colnames(env_vals) <- names(env)

cat("NAs per variable before filtering:\n")
print(colSums(is.na(env_vals)))

# Step 1: Pearson correlation filtering (threshold = 0.7)
# One variable from each pair with |r| > 0.7 is removed.
cor_mat <- cor(env_vals, use = "complete.obs")

corrplot(cor_mat, method = "color", tl.cex = 0.7,
         title = "Pearson correlation — 19 bioclimatic variables",
         mar   = c(0, 0, 2, 0))

high_cor  <- findCorrelation(cor_mat, cutoff = 0.7)
env_uncor <- env[[-high_cor]]

cat("Variables retained after correlation filtering:", nlyr(env_uncor), "\n")
cat("Retained:", paste(names(env_uncor), collapse = ", "), "\n")

# Step 2: Variance Inflation Factor filtering (threshold = 10)
# Removes variables with VIF >= 10 iteratively.
vif_res <- vifstep(env_uncor)
env_vif <- exclude(env_uncor, vif_res)

cat("Variables retained after VIF filtering:", nlyr(env_vif), "\n")
cat("Retained:", paste(names(env_vif), collapse = ", "), "\n")
cat("\nVIF values:\n")
print(vif_res)

# Step 3: Retain all statistically filtered variables
# All 8 variables passed both filters (max VIF = 3.90, all |r| < 0.7)
# representing non-redundant thermal and precipitation dimensions:
#
# bio5  = Max Temperature of Warmest Month   — upper thermal ceiling
# bio9  = Mean Temperature of Driest Quarter — thermal stress during drought
# bio10 = Mean Temperature of Warmest Quarter — warm season thermal environment
# bio11 = Mean Temperature of Coldest Quarter — cold tolerance / invasion barrier
# bio12 = Annual Precipitation               — overall water availability
# bio13 = Precipitation of Wettest Month     — peak moisture availability
# bio18 = Precipitation of Warmest Quarter   — moisture during thermal stress
# bio19 = Precipitation of Coldest Quarter   — winter moisture regime

final_vars <- c("bio5", "bio9", "bio10", "bio11",
                "bio12", "bio13", "bio18", "bio19")
env_final  <- env[[final_vars]]

cat("Final variable set:", paste(names(env_final), collapse = ", "), "\n")

# =========================================
# SECTION 7: NA CHECK AND CORRECTION
# =========================================
# NA values in WorldClim 2.1 at occurrence points arise from
# land-mask discontinuities at coastal grid cells.
# Where detected, NAs are filled using a 3x3 focal mean
# restricted to NA cells only (standard SDM practice).

na_check <- extract(env_final,
                    as.matrix(occ_thin[, c("lon", "lat")]))[, names(env_final)]
colnames(na_check) <- final_vars

cat("NA count per variable at thinned occurrence localities:\n")
print(colSums(is.na(na_check)))

vars_with_na <- names(which(colSums(is.na(na_check)) > 0))

if (length(vars_with_na) == 0) {
  cat("No NAs at occurrence localities. Checking background compatibility.\n")
} else {
  cat("Applying focal fill to:", paste(vars_with_na, collapse = ", "), "\n")
  for (v in vars_with_na) {
    filled         <- focal(env_final[[v]], w = 3, fun = "mean",
                            na.policy = "only", na.rm = TRUE)
    env_final[[v]] <- filled
    cat("Focal fill applied to", v, "— remaining NAs:",
        global(is.na(env_final[[v]]), sum)[[1]], "\n")
  }
}

# =========================================
# SECTION 8: GLOBAL MODEL (ALL OCCURRENCES)
# =========================================

# Generate background points from bio5 land mask
# (ensures all background points fall on land cells with valid data)
bg <- spatSample(env_final[["bio5"]], size = 10000, method = "random",
                 na.rm = TRUE, xy = TRUE)
bg <- as.data.frame(bg)[, 1:2]
colnames(bg) <- c("lon", "lat")

cat("Background points generated:", nrow(bg), "\n")

# Combine presences and background
pres    <- occ_thin
pres$pa <- 1
bg$pa   <- 0

dat <- rbind(pres[, c("lon", "lat", "pa")], bg)

# Extract environmental values by name (not position) to prevent silent drops
coords   <- as.matrix(dat[, c("lon", "lat")])
env_vals <- extract(env_final, coords)[, names(env_final)]
env_vals <- as.data.frame(env_vals)

# Final combined dataset
dat_full <- na.omit(cbind(dat, env_vals))

cat("Final dataset — presences:", sum(dat_full$pa == 1), "\n")
cat("Final dataset — background:", sum(dat_full$pa == 0), "\n")
cat("Variables in final model:",
    paste(names(dat_full)[4:ncol(dat_full)], collapse = ", "), "\n")

# Fit MaxNet model (linear + quadratic + hinge features)
model_all <- maxnet(
  p       = dat_full$pa,
  data    = dat_full[, final_vars],
  f       = maxnet.formula(dat_full$pa, dat_full[, final_vars], classes = "lqh"),
  regmult = 0.5
)

# Predict global suitability
pred_global <- predict(env_final, model_all,
                       fun   = function(m, d) predict(m, d, type = "cloglog"),
                       na.rm = TRUE)

# Plot global prediction
plot(pred_global,
     main = expression(italic("Zaprionus indianus") ~ "— Global Habitat Suitability"),
     col  = rev(terrain.colors(100)))

# Save global prediction raster
writeRaster(pred_global,
            "C:/Users/Avirup/OneDrive/Desktop/Zap_SDM/outputs/global_suitability.tif",
            overwrite = TRUE)

# Predict for India
india      <- geodata::gadm("IND", level = 0, path = tempdir())
env_india  <- mask(crop(env_final, india), india)
pred_india <- predict(env_india, model_all,
                      fun   = function(m, d) predict(m, d, type = "cloglog"),
                      na.rm = TRUE)

# Plot India prediction
plot(pred_india,
     main = expression(italic("Zaprionus indianus") ~ "— India Habitat Suitability"),
     col  = rev(terrain.colors(100)))
points(occ_thin, col = "red", pch = 20, cex = 0.5)

# Save India prediction raster
writeRaster(pred_india,
            "C:/Users/Avirup/OneDrive/Desktop/Zap_SDM/outputs/india_suitability.tif",
            overwrite = TRUE)

# =========================================
# SECTION 9: NATIVE RANGE MODEL
# =========================================
# Native range defined as the Afrotropical region
# Longitude: -20 to 55, Latitude: -35 to 38

occ_thin$region <- ifelse(
  occ_thin$lon >= -20 & occ_thin$lon <= 55 &
    occ_thin$lat >= -35 & occ_thin$lat <= 38,
  "native", "invasive"
)

occ_native   <- subset(occ_thin, region == "native")
occ_invasive <- subset(occ_thin, region == "invasive")

cat("Native occurrences:", nrow(occ_native), "\n")
cat("Invasive occurrences:", nrow(occ_invasive), "\n")

# Build native-only model
pres_native    <- occ_native[, c("lon", "lat")]
pres_native$pa <- 1
bg_use         <- bg[, c("lon", "lat")]
bg_use$pa      <- 0

dat_nat  <- rbind(pres_native, bg_use)
coords   <- as.matrix(dat_nat[, c("lon", "lat")])
env_nat  <- as.data.frame(extract(env_final, coords)[, names(env_final)])

dat_nat_full <- na.omit(cbind(dat_nat, env_nat))

X <- data.frame(lapply(dat_nat_full[, final_vars], as.numeric))
y <- dat_nat_full$pa

model_nat <- maxnet(
  p       = y,
  data    = X,
  f       = maxnet.formula(y, X, classes = "lqh"),
  regmult = 0.5
)

# Predict native model globally and for India
pred_nat_global <- predict(env_final, model_nat,
                           fun   = function(m, d) predict(m, d, type = "cloglog"),
                           na.rm = TRUE)
plot(pred_nat_global,
     main = expression(italic("Zaprionus indianus") ~ "— Global Habitat Suitability (Native)"),
     col  = rev(terrain.colors(100)))

pred_nat_india <- predict(env_india, model_nat,
                          fun   = function(m, d) predict(m, d, type = "cloglog"),
                          na.rm = TRUE)

plot(pred_nat_india,
     main = expression(italic("Zaprionus indianus") ~ "— Global Habitat Suitability (Native)"),
     col  = rev(terrain.colors(100)))

# Save native model predictions
writeRaster(pred_nat_global,
            "C:/Users/Avirup/OneDrive/Desktop/Zap_SDM/outputs/native_global_suitability.tif",
            overwrite = TRUE)

writeRaster(pred_nat_india,
            "C:/Users/Avirup/OneDrive/Desktop/Zap_SDM/outputs/native_india_suitability.tif",
            overwrite = TRUE)

# =========================================
# SECTION 10: RESPONSE CURVES
# =========================================

# Ensure pROC functions take precedence over PresenceAbsence
conflictRules("pROC", prefer = c("roc", "auc", "coords"))

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

for (v in final_vars) {

  x <- seq(min(dat_full[[v]], na.rm = TRUE),
           max(dat_full[[v]], na.rm = TRUE),
           length.out = 100)

  newdata           <- as.data.frame(
                         matrix(rep(colMeans(dat_full[, final_vars], na.rm = TRUE),
                                    each = 100), nrow = 100))
  colnames(newdata) <- final_vars
  newdata[[v]]      <- x

  y_full <- predict(model_all, newdata, type = "cloglog")
  y_nat  <- predict(model_nat, newdata,  type = "cloglog")

  plot(x, y_full,
       type = "l", lwd = 2, col = "red",
       ylim = c(0, max(y_full, y_nat, na.rm = TRUE)),
       xlab = v, ylab = "Suitability", main = v)

  lines(x, y_nat, col = "blue", lwd = 2)
}

# Legend on final panel
plot.new()
legend("center",
       legend = c("Global model", "Native model"),
       col    = c("red", "blue"),
       lwd    = 2, bty = "n", cex = 1.1)

# =========================================
# SECTION 11: MODEL EVALUATION (AUC + TSS)
# =========================================

# --- 11a: Full training data evaluation ---
pred_vals <- as.numeric(predict(model_all, dat_full[, final_vars],
                                type = "cloglog"))

roc_obj  <- pROC::roc(dat_full$pa, pred_vals, quiet = TRUE)
auc_full <- as.numeric(pROC::auc(roc_obj))

opt       <- pROC::coords(roc_obj, x = "best",
                          ret = c("threshold", "sensitivity", "specificity"))
threshold <- as.numeric(opt["threshold"])
pred_bin  <- ifelse(pred_vals >= threshold, 1, 0)

TP <- sum(dat_full$pa == 1 & pred_bin == 1)
TN <- sum(dat_full$pa == 0 & pred_bin == 0)
FP <- sum(dat_full$pa == 0 & pred_bin == 1)
FN <- sum(dat_full$pa == 1 & pred_bin == 0)

TSS_full <- (TP / (TP + FN)) + (TN / (TN + FP)) - 1

cat("Training AUC:", round(auc_full, 3), "\n")
cat("Training TSS:", round(TSS_full, 3), "\n")

# --- 11b: Held-out test data evaluation (70/30 split) ---
set.seed(42)

train_idx  <- sample(1:nrow(dat_full), size = 0.7 * nrow(dat_full))
train_data <- dat_full[train_idx, ]
test_data  <- dat_full[-train_idx, ]

model_eval <- maxnet(
  p       = train_data$pa,
  data    = train_data[, final_vars],
  f       = maxnet.formula(train_data$pa, train_data[, final_vars],
                           classes = "lqh"),
  regmult = 0.5
)

pred_test <- as.numeric(predict(model_eval, test_data[, final_vars],
                                type = "cloglog"))

roc_test   <- pROC::roc(test_data$pa, pred_test, quiet = TRUE)
auc_test   <- as.numeric(pROC::auc(roc_test))

opt_test    <- pROC::coords(roc_test, x = "best",
                            ret = c("threshold", "sensitivity", "specificity"))
thresh_test <- as.numeric(opt_test["threshold"])
pred_bin_t  <- ifelse(pred_test >= thresh_test, 1, 0)

TP_t <- sum(test_data$pa == 1 & pred_bin_t == 1)
TN_t <- sum(test_data$pa == 0 & pred_bin_t == 0)
FP_t <- sum(test_data$pa == 0 & pred_bin_t == 1)
FN_t <- sum(test_data$pa == 1 & pred_bin_t == 0)

TSS_test <- (TP_t / (TP_t + FN_t)) + (TN_t / (TN_t + FP_t)) - 1

cat("Test AUC (70/30 split):", round(auc_test, 3), "\n")
cat("Test TSS (70/30 split):", round(TSS_test, 3), "\n")

# =========================================
# SECTION 12: BINARY SUITABILITY MAPS
# =========================================
# MaxTSS threshold derived from full model predictions
# at occurrence and background localities

pred_pres <- as.numeric(extract(pred_global,
                                occ_thin[, c("lon", "lat")])[, 2])
pred_bg   <- as.numeric(extract(pred_global,
                                bg[, c("lon", "lat")])[, 2])

obs_all  <- c(rep(1, length(pred_pres)), rep(0, length(pred_bg)))
pred_all <- c(pred_pres, pred_bg)

valid    <- !is.na(pred_all)
obs_all  <- obs_all[valid]
pred_all <- pred_all[valid]

roc_full        <- roc(obs_all, pred_all, quiet = TRUE)
opt_full        <- coords(roc_full, x = "best",
                          ret = c("threshold", "sensitivity", "specificity"))
threshold_final <- as.numeric(opt_full["threshold"])

cat("MaxTSS threshold for binary maps:", round(threshold_final, 3), "\n")

binary_global <- pred_global >= threshold_final
binary_india  <- pred_india  >= threshold_final

# Save binary rasters
writeRaster(binary_global,
            "C:/Users/Avirup/OneDrive/Desktop/Zap_SDM/outputs/binary_global.tif",
            overwrite = TRUE)
writeRaster(binary_india,
            "C:/Users/Avirup/OneDrive/Desktop/Zap_SDM/outputs/binary_india.tif",
            overwrite = TRUE)

# --- Plot settings ---
col_unsuitable <- "#D9D9D9"
col_suitable   <- "#2166AC"
world          <- geodata::world(path = tempdir())

# Global binary map
par(mfrow = c(1, 1), mar = c(3, 3, 3, 3))

plot(binary_global,
     col      = c(col_unsuitable, col_suitable),
     main     = expression(italic("Zaprionus indianus") ~
                           "— Global Habitat Suitability"),
     cex.main = 1.1,
     legend   = FALSE,
     axes     = TRUE,
     box      = FALSE)
lines(world, col = "white", lwd = 0.4)
legend("bottomleft",
       legend  = c("Unsuitable", "Suitable"),
       fill    = c(col_unsuitable, col_suitable),
       border  = "grey40",
       title   = expression(bold("Suitability")),
       bty     = "o",
       box.col = "grey40",
       bg      = "white",
       cex     = 0.85,
       inset   = c(0.01, 0.02))

# India binary map
par(mfrow = c(1, 1), mar = c(3, 3, 3, 3))

plot(binary_india,
     col      = c(col_unsuitable, col_suitable),
     main     = expression(italic("Zaprionus indianus") ~
                           "— Habitat Suitability in India"),
     cex.main = 1.1,
     legend   = FALSE,
     axes     = TRUE,
     box      = FALSE)
lines(india, col = "grey30", lwd = 0.7)
legend("bottomleft",
       legend  = c("Unsuitable", "Suitable"),
       fill    = c(col_unsuitable, col_suitable),
       border  = "grey40",
       title   = expression(bold("Suitability")),
       bty     = "o",
       box.col = "grey40",
       bg      = "white",
       cex     = 0.85,
       inset   = c(0.01, 0.02))

# =========================================
# SECTION 13: CLIMATIC NICHE ANALYSIS
# =========================================
# Native range = Africa (Afrotropical region)
# Invaded range = rest of the world
# Regional backgrounds sampled separately to represent
# available climatic space in each region

africa_countries <- c(
  "Algeria", "Angola", "Benin", "Botswana", "Burkina Faso", "Burundi",
  "Cameroon", "Cape Verde", "Central African Republic", "Chad", "Comoros",
  "Congo", "Democratic Republic of the Congo", "Djibouti", "Egypt",
  "Equatorial Guinea", "Eritrea", "Eswatini", "Ethiopia", "Gabon", "Gambia",
  "Ghana", "Guinea", "Guinea-Bissau", "Ivory Coast", "Kenya", "Lesotho",
  "Liberia", "Libya", "Madagascar", "Malawi", "Mali", "Mauritania", "Mauritius",
  "Morocco", "Mozambique", "Namibia", "Niger", "Nigeria", "Rwanda",
  "Sao Tome and Principe", "Senegal", "Seychelles", "Sierra Leone", "Somalia",
  "South Africa", "South Sudan", "Sudan", "Tanzania", "Togo", "Tunisia",
  "Uganda", "Zambia", "Zimbabwe"
)

africa     <- world[world$NAME_0 %in% africa_countries, ]
rest_world <- world[!(world$NAME_0 %in% africa_countries), ]

# Sample regional backgrounds
env_nat_region <- terra::mask(env_final, africa)
env_inv_region <- terra::mask(env_final, rest_world)

bg_nat <- as.data.frame(terra::spatSample(env_nat_region, size = 10000,
                                           na.rm = TRUE, xy = TRUE))
colnames(bg_nat)[1:2] <- c("lon", "lat")
bg_nat$species <- 0

bg_inv <- as.data.frame(terra::spatSample(env_inv_region, size = 10000,
                                           na.rm = TRUE, xy = TRUE))
colnames(bg_inv)[1:2] <- c("lon", "lat")
bg_inv$species <- 0

# Build presence dataframes for each region
nat_pres <- occ_native[, c("lon", "lat")]
nat_pres <- cbind(nat_pres,
                  extract(env_final,
                          as.matrix(nat_pres))[, names(env_final)])
nat_pres$species <- 1

inv_pres <- occ_invasive[, c("lon", "lat")]
inv_pres <- cbind(inv_pres,
                  extract(env_final,
                          as.matrix(inv_pres))[, names(env_final)])
inv_pres$species <- 1

# Extract env values for background points
bg_nat <- cbind(bg_nat[, c("lon", "lat")],
                extract(env_final,
                        as.matrix(bg_nat[, c("lon", "lat")]))[, names(env_final)])
bg_nat$species <- 0

bg_inv <- cbind(bg_inv[, c("lon", "lat")],
                extract(env_final,
                        as.matrix(bg_inv[, c("lon", "lat")]))[, names(env_final)])
bg_inv$species <- 0

# Combine and save
native_df   <- na.omit(rbind(nat_pres, bg_nat))
invasive_df <- na.omit(rbind(inv_pres, bg_inv))

write.csv(native_df,   "C:/Users/Avirup/OneDrive/Desktop/Zap_SDM/outputs/ecospat_native.csv",
          row.names = FALSE)
write.csv(invasive_df, "C:/Users/Avirup/OneDrive/Desktop/Zap_SDM/outputs/ecospat_invasive.csv",
          row.names = FALSE)

cat("Native ecospat dataset rows:", nrow(native_df), "\n")
cat("Invasive ecospat dataset rows:", nrow(invasive_df), "\n")

# =========================================
# SECTION 14: ECOSPAT NICHE ANALYSIS
# =========================================

nat <- read.csv("C:/Users/Avirup/OneDrive/Desktop/Zap_SDM/outputs/ecospat_native.csv")
inv <- read.csv("C:/Users/Avirup/OneDrive/Desktop/Zap_SDM/outputs/ecospat_invasive.csv")

# Environmental columns (all bio variables, excluding lon, lat, species)
env_cols <- names(nat)[!names(nat) %in% c("lon", "lat", "species")]

# PCA on pooled environmental space (native + invaded backgrounds)
pca.env <- dudi.pca(rbind(nat, inv)[, env_cols],
                    scannf = FALSE, nf = 2)

ecospat.plot.contrib(contrib = pca.env$co, eigen = pca.env$eig)

# Project scores onto PCA axes
scores.globclim <- pca.env$li

scores.sp.nat   <- suprow(pca.env,
                           nat[nat$species == 1, env_cols])$li
scores.sp.inv   <- suprow(pca.env,
                           inv[inv$species == 1, env_cols])$li
scores.clim.nat <- suprow(pca.env, nat[, env_cols])$li
scores.clim.inv <- suprow(pca.env, inv[, env_cols])$li

# Build dynamic occurrence grids
grid.clim.nat <- ecospat.grid.clim.dyn(glob  = scores.globclim,
                                        glob1 = scores.clim.nat,
                                        sp    = scores.sp.nat,
                                        R     = 100, th.sp = 0)

grid.clim.inv <- ecospat.grid.clim.dyn(glob  = scores.globclim,
                                        glob1 = scores.clim.inv,
                                        sp    = scores.sp.inv,
                                        R     = 100, th.sp = 0)

# Schoener's D niche overlap
D.overlap <- ecospat.niche.overlap(grid.clim.nat, grid.clim.inv,
                                   cor = TRUE)$D
cat("Schoener's D:", round(D.overlap, 3), "\n")

# Niche equivalency test (100 permutations)
eq.test <- ecospat.niche.equivalency.test(grid.clim.nat, grid.clim.inv,
                                           rep = 100)
ecospat.plot.overlap.test(eq.test, "D", "Equivalency")

# Niche similarity test (1000 permutations)
sim.test <- ecospat.niche.similarity.test(grid.clim.nat, grid.clim.inv,
                                           rep = 1000, rand.type = 2)
ecospat.plot.overlap.test(sim.test, "D", "Similarity")

# Niche dynamics: stability, expansion, unfilling
niche.dyn <- ecospat.niche.dyn.index(grid.clim.nat, grid.clim.inv)
cat("\nNiche dynamics:\n")
print(niche.dyn)

# Plot niche dynamics
ecospat.plot.niche.dyn(grid.clim.nat, grid.clim.inv,
                        quant      = 0.25,
                        interest   = 2,
                        title      = "",
                        name.axis1 = "PC1",
                        name.axis2 = "PC2")

# Centroid shift
ecospat.shift.centroids(scores.sp.nat, scores.sp.inv,
                         scores.clim.nat, scores.clim.inv,
                         col = "black")

# =========================================
# SECTION 15: UNIVARIATE NICHE PLOTS
# =========================================

nat_pres_plot <- nat[nat$species == 1, ]
inv_pres_plot <- inv[inv$species == 1, ]

plot_var <- function(var) {
  d_nat_bg <- density(nat[[var]],          na.rm = TRUE)
  d_inv_bg <- density(inv[[var]],          na.rm = TRUE)
  d_nat    <- density(nat_pres_plot[[var]], na.rm = TRUE)
  d_inv    <- density(inv_pres_plot[[var]], na.rm = TRUE)

  plot(d_nat_bg,
       col  = "grey40", lwd = 2,
       main = var, xlab = var, ylab = "Density",
       ylim = c(0, max(d_nat_bg$y, d_inv_bg$y,
                       d_nat$y,    d_inv$y, na.rm = TRUE)))

  lines(d_inv_bg, col = "grey40", lwd = 2, lty = 2)
  lines(d_nat,    col = "#2ecc71", lwd = 2)
  lines(d_inv,    col = "#e74c3c", lwd = 2)

  x_common <- seq(min(d_nat$x, d_inv$x),
                  max(d_nat$x, d_inv$x), length = 500)
  y_nat_i  <- approx(d_nat$x, d_nat$y, xout = x_common)$y
  y_inv_i  <- approx(d_inv$x, d_inv$y, xout = x_common)$y

  polygon(c(x_common, rev(x_common)),
          c(pmin(y_nat_i, y_inv_i), rep(0, length(x_common))),
          col    = adjustcolor("purple", 0.4),
          border = NA)
}

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

for (v in env_cols) {
  plot_var(v)
}

# Legend on final panel
plot.new()
legend("center",
       legend = c("Native climate space", "Invasive climate space",
                  "Native niche", "Invasive niche", "Overlap"),
       col    = c("grey40", "grey40", "#2ecc71", "#e74c3c", "purple"),
       lwd    = 2,
       lty    = c(1, 2, 1, 1, 1),
       bty    = "n", cex = 0.9)

# =========================================
# SECTION 16: MESS ANALYSIS
# =========================================
# MESS (Multivariate Environmental Similarity Surface) identifies
# regions in the invaded range that fall outside the environmental
# space of the native range (novel environments / extrapolation risk).
# Negative MESS values indicate novel climatic conditions.

# Build MESS input from native and invasive ecospat datasets
# Calibration = native range occurrences + background
# Projection  = invasive range occurrences + background

mess_cols <- c("lon", "lat", env_cols)

cal  <- na.omit(native_df[, mess_cols])
proj <- na.omit(invasive_df[, mess_cols])

cat("MESS calibration points (native):", nrow(cal), "\n")
cat("MESS projection points (invasive):", nrow(proj), "\n")

# Compute MESS
mess.object <- ecospat.mess(proj, cal, w = "default")

# Base plot
par(mfrow = c(1, 3))
ecospat.plot.mess(mess.object = mess.object)

# ggplot2 MESS maps
mess.df  <- as.data.frame(mess.object)
colname5 <- colnames(mess.df)[5]

p1 <- ggplot(mess.df, aes(x = lon, y = lat, fill = MESS)) +
  geom_tile() +
  scale_fill_gradient2(low      = "red",
                       mid      = "white",
                       high     = "blue",
                       midpoint = 0,
                       name     = "MESS") +
  labs(title = "MESS — multivariate similarity",
       x = "Longitude", y = "Latitude") +
  theme_minimal() +
  theme(legend.position = "right")

p2 <- ggplot(mess.df, aes(x = lon, y = lat, fill = MESSw)) +
  geom_tile() +
  scale_fill_gradient2(low      = "red",
                       mid      = "white",
                       high     = "blue",
                       midpoint = 0,
                       name     = "MESSw") +
  labs(title = "MESSw — weighted similarity",
       x = "Longitude", y = "Latitude") +
  theme_minimal() +
  theme(legend.position = "right")

p3 <- ggplot(mess.df, aes(x = lon, y = lat,
                           fill = .data[[colname5]])) +
  geom_tile() +
  scale_fill_gradient(low    = "white",
                      high   = "red",
                      limits = c(0, length(env_cols)),
                      name   = "# Novel\nVariables") +
  labs(title = "# Variables with novel conditions",
       x = "Longitude", y = "Latitude") +
  theme_minimal() +
  theme(legend.position = "right")

grid.arrange(p1, p2, p3, ncol = 2)

# Save MESS plots
ggsave("C:/Users/Avirup/OneDrive/Desktop/Zap_SDM/outputs/MESS_plots.pdf",
       plot     = grid.arrange(p1, p2, p3, ncol = 2),
       width    = 12,
       height   = 8,
       dpi      = 300)

cat("\nPipeline complete. All outputs saved to:\n")
cat("C:/Users/Avirup/OneDrive/Desktop/Zap_SDM/outputs/\n")

mess.df  <- as.data.frame(mess.object)
colname5 <- colnames(mess.df)[5]

# Get world boundaries for context
world <- geodata::world(path = tempdir())
world_df <- as.data.frame(world, geom = "XY")

p1 <- ggplot() +
  geom_point(data = mess.df,
             aes(x = lon, y = lat, color = MESS),
             size = 0.8, alpha = 0.7) +
  scale_color_gradient2(low      = "#d73027",
                        mid      = "white",
                        high     = "#4575b4",
                        midpoint = 0,
                        name     = "MESS",
                        limits   = c(-231, 62),
                        oob      = scales::squish) +
  geom_sf(data = sf::st_as_sf(world),
          fill = NA, color = "grey40", linewidth = 0.2) +
  coord_sf(xlim = c(-180, 180), ylim = c(-60, 80)) +
  labs(title    = "Multivariate Environmental Similarity Surface (MESS)",
       subtitle = paste0("Blue = similar to native range  |  ",
                         "Red = novel conditions  |  ",
                         "54.6% of invaded range is novel (MESS < 0)"),
       x = NULL, y = NULL) +
  theme_minimal(base_size = 11) +
  theme(legend.position  = "right",
        panel.grid.major = element_line(color = "grey90", linewidth = 0.3),
        plot.title       = element_text(face = "bold", size = 12),
        plot.subtitle    = element_text(size = 9, color = "grey40"))

p2 <- ggplot() +
  geom_point(data = mess.df,
             aes(x = lon, y = lat, color = MESSw),
             size = 0.8, alpha = 0.7) +
  scale_color_gradient2(low      = "#d73027",
                        mid      = "white",
                        high     = "#4575b4",
                        midpoint = 0,
                        name     = "MESSw") +
  geom_sf(data = sf::st_as_sf(world),
          fill = NA, color = "grey40", linewidth = 0.2) +
  coord_sf(xlim = c(-180, 180), ylim = c(-60, 80)) +
  labs(title = "Weighted MESS (MESSw)",
       x = NULL, y = NULL) +
  theme_minimal(base_size = 11) +
  theme(legend.position  = "right",
        panel.grid.major = element_line(color = "grey90", linewidth = 0.3),
        plot.title       = element_text(face = "bold", size = 12))

p3 <- ggplot() +
  geom_point(data = mess.df,
             aes(x = lon, y = lat,
                 color = .data[[colname5]]),
             size = 0.8, alpha = 0.7) +
  scale_color_gradient(low   = "#ffffcc",
                       high  = "#d73027",
                       name  = "# Novel\nVariables",
                       limits = c(0, 8),
                       breaks = c(0, 2, 4, 6, 8)) +
  geom_sf(data = sf::st_as_sf(world),
          fill = NA, color = "grey40", linewidth = 0.2) +
  coord_sf(xlim = c(-180, 180), ylim = c(-60, 80)) +
  labs(title    = "Number of variables with novel conditions",
       subtitle = "Higher values indicate greater extrapolation risk",
       x = NULL, y = NULL) +
  theme_minimal(base_size = 11) +
  theme(legend.position  = "right",
        panel.grid.major = element_line(color = "grey90", linewidth = 0.3),
        plot.title       = element_text(face = "bold", size = 12),
        plot.subtitle    = element_text(size = 9, color = "grey40"))

# Combine
mess_combined <- grid.arrange(p1, p2, p3, ncol = 1,
                              top = grid::textGrob(
                                expression(italic("Zaprionus indianus") ~
                                             "— Climatic Novelty in the Invaded Range"),
                                gp = grid::gpar(fontsize = 14, fontface = "bold")))

# Save
ggsave("C:/Users/Avirup/OneDrive/Desktop/Zap_SDM/outputs/MESS_final.pdf",
       plot   = mess_combined,
       width  = 12,
       height = 14,
       dpi    = 300)

ggsave("C:/Users/Avirup/OneDrive/Desktop/Zap_SDM/outputs/MESS_final.png",
       plot   = mess_combined,
       width  = 12,
       height = 14,
       dpi    = 300,
       bg     = "white")

cat("MESS figures saved.\n")

# Add Africa outline to show native range reference
africa_sf <- sf::st_as_sf(africa)

p1 <- p1 +
  geom_sf(data   = africa_sf,
          fill   = "grey85",
          color  = "grey40",
          linewidth = 0.3,
          alpha  = 0.5) +
  annotate("text", x = 20, y = 5, label = "Native range",
           size = 2.8, color = "grey30", fontface = "italic")

p2 <- p2 +
  geom_sf(data   = africa_sf,
          fill   = "grey85",
          color  = "grey40",
          linewidth = 0.3,
          alpha  = 0.5) +
  annotate("text", x = 20, y = 5, label = "Native range",
           size = 2.8, color = "grey30", fontface = "italic")

p3 <- p3 +
  geom_sf(data   = africa_sf,
          fill   = "grey85",
          color  = "grey40",
          linewidth = 0.3,
          alpha  = 0.5) +
  annotate("text", x = 20, y = 5, label = "Native range",
           size = 2.8, color = "grey30", fontface = "italic")

# Replot
mess_combined <- grid.arrange(p1, p2, p3, ncol = 1,
                              top = grid::textGrob(
                                expression(italic("Zaprionus indianus") ~
                                             "— Climatic Novelty in the Invaded Range"),
                                gp = grid::gpar(fontsize = 14, fontface = "bold")))

ggsave("C:/Users/Avirup/OneDrive/Desktop/Zap_SDM/outputs/MESS_final.png",
       plot   = mess_combined,
       width  = 12,
       height = 14,
       dpi    = 300,
       bg     = "white")



# =========================================
# FIGURE CODE — MAIN TEXT AND SUPPLEMENTARY
# Zaprionus indianus SDM
# =========================================

library(sf)
library(ggplot2)
library(gridExtra)
library(terra)
library(RColorBrewer)
library(scales)

# Output directories
dir.create("C:/Users/Avirup/OneDrive/Desktop/Zap_SDM/outputs/figures/main",
           showWarnings = FALSE, recursive = TRUE)
dir.create("C:/Users/Avirup/OneDrive/Desktop/Zap_SDM/outputs/figures/supplementary",
           showWarnings = FALSE, recursive = TRUE)

# Shared objects
world_sf  <- sf::st_as_sf(geodata::world(path = tempdir()))
india_sf  <- sf::st_as_sf(geodata::gadm("IND", level = 1, path = tempdir()))
africa_sf <- sf::st_as_sf(africa)

# Colour palettes
pal_suit    <- colorRampPalette(c("#f7f7f7", "#d1e5f0", "#4393c3",
                                  "#2166ac", "#053061"))(100)
pal_mess    <- colorRampPalette(c("#d73027", "white", "#4575b4"))(100)
col_native  <- "#2ecc71"
col_invasive <- "#e74c3c"

# =========================================
# MAIN TEXT FIGURE 1
# Global occurrence map — native vs invasive
# =========================================

occ_plot <- rbind(
  data.frame(occ_native[, c("lon","lat")],   range = "Native"),
  data.frame(occ_invasive[, c("lon","lat")], range = "Invasive")
)

fig1 <- ggplot() +
  geom_sf(data = world_sf, fill = "grey92", color = "grey60",
          linewidth = 0.2) +
  geom_sf(data = africa_sf, fill = "grey80", color = "grey50",
          linewidth = 0.25) +
  geom_point(data = occ_plot,
             aes(x = lon, y = lat, color = range, shape = range),
             size = 1.8, alpha = 0.8) +
  scale_color_manual(values = c("Native" = col_native,
                                "Invasive" = col_invasive),
                     name = "Range") +
  scale_shape_manual(values = c("Native" = 17, "Invasive" = 16),
                     name = "Range") +
  coord_sf(xlim = c(-180, 180), ylim = c(-90, 90)) +
  labs(title    = expression(italic("Zaprionus indianus") ~
                               "— Global occurrence records"),
       subtitle = paste0("Native: n = ", nrow(occ_native),
                         "  |  Invasive: n = ", nrow(occ_invasive),
                         "  |  Total thinned: n = ", nrow(occ_thin)),
       x = NULL, y = NULL) +
  theme_minimal(base_size = 11) +
  theme(legend.position   = "bottom",
        legend.title      = element_text(face = "bold"),
        panel.grid.major  = element_line(color = "grey88", linewidth = 0.3),
        plot.title        = element_text(face = "bold", size = 13),
        plot.subtitle     = element_text(size = 9, color = "grey40"))

ggsave("C:/Users/Avirup/OneDrive/Desktop/Zap_SDM/outputs/figures/main/Fig1_occurrences.pdf",
       fig1, width = 10, height = 6, dpi = 300)
ggsave("C:/Users/Avirup/OneDrive/Desktop/Zap_SDM/outputs/figures/main/Fig1_occurrences.png",
       fig1, width = 10, height = 6, dpi = 300, bg = "white")

cat("Figure 1 saved.\n")

# =========================================
# MAIN TEXT FIGURE 2
# Global habitat suitability — continuous
# Panel A: Global model | Panel B: Native model
# =========================================

# Convert rasters to dataframes
pred_global_df <- as.data.frame(pred_global, xy = TRUE, na.rm = TRUE)
colnames(pred_global_df) <- c("lon", "lat", "suitability")

pred_nat_df <- as.data.frame(pred_nat_global, xy = TRUE, na.rm = TRUE)
colnames(pred_nat_df) <- c("lon", "lat", "suitability")

make_global_suit_map <- function(df, title_text, subtitle_text) {
  ggplot() +
    geom_tile(data = df, aes(x = lon, y = lat, fill = suitability)) +
    scale_fill_gradientn(colors = pal_suit,
                         limits = c(0, 1),
                         name   = "Suitability",
                         breaks = c(0, 0.25, 0.5, 0.75, 1.0),
                         labels = c("0", "0.25", "0.50", "0.75", "1.0")) +
    geom_sf(data = world_sf, fill = NA, color = "grey50", linewidth = 0.2) +
    geom_point(data = occ_thin, aes(x = lon, y = lat),
               color = "white", size = 0.4, alpha = 0.6, shape = 3) +
    coord_sf(xlim = c(-180, 180), ylim = c(-90, 90)) +
    labs(title = title_text, subtitle = subtitle_text,
         x = NULL, y = NULL) +
    theme_minimal(base_size = 11) +
    theme(legend.position  = "right",
          panel.grid.major = element_line(color = "grey88", linewidth = 0.3),
          plot.title       = element_text(face = "bold", size = 12),
          plot.subtitle    = element_text(size = 9, color = "grey40"))
}

fig2a <- make_global_suit_map(
  pred_global_df,
  expression(bold("A") ~ " Global model — all occurrences"),
  "Trained on native + invasive range occurrences"
)

fig2b <- make_global_suit_map(
  pred_nat_df,
  expression(bold("B") ~ " Native range model"),
  "Trained on Afrotropical occurrences only"
)

fig2 <- grid.arrange(fig2a, fig2b, ncol = 1,
                     top = grid::textGrob(
                       expression(italic("Zaprionus indianus") ~
                                    "— Global Habitat Suitability"),
                       gp = grid::gpar(fontsize = 14, fontface = "bold")))

ggsave("C:/Users/Avirup/OneDrive/Desktop/Zap_SDM/outputs/figures/main/Fig2_global_suitability.pdf",
       fig2, width = 10, height = 10, dpi = 300)
ggsave("C:/Users/Avirup/OneDrive/Desktop/Zap_SDM/outputs/figures/main/Fig2_global_suitability.png",
       fig2, width = 10, height = 10, dpi = 300, bg = "white")

cat("Figure 2 saved.\n")

# =========================================
# MAIN TEXT FIGURE 3
# India habitat suitability
# Panel A: Global model | Panel B: Native model
# =========================================

pred_india_df <- as.data.frame(pred_india, xy = TRUE, na.rm = TRUE)
colnames(pred_india_df) <- c("lon", "lat", "suitability")

pred_nat_india_df <- as.data.frame(pred_nat_india, xy = TRUE, na.rm = TRUE)
colnames(pred_nat_india_df) <- c("lon", "lat", "suitability")

occ_india_pts <- occ_thin[occ_thin$lon >= 68 & occ_thin$lon <= 97 &
                            occ_thin$lat >= 8  & occ_thin$lat <= 37, ]

make_india_map <- function(df, title_text, subtitle_text) {
  ggplot() +
    geom_raster(data = df, aes(x = lon, y = lat, fill = suitability)) +
    scale_fill_gradientn(colors = pal_suit,
                         limits = c(0, 1),
                         name   = "Suitability",
                         breaks = c(0, 0.25, 0.5, 0.75, 1.0),
                         labels = c("0", "0.25", "0.50", "0.75", "1.0")) +
    geom_sf(data = india_sf, fill = NA, color = "grey40", linewidth = 0.3) +
    geom_point(data = occ_india_pts, aes(x = lon, y = lat),
               color = "white", size = 1.5, shape = 3, stroke = 0.8) +
    coord_sf(xlim = c(68, 97), ylim = c(8, 37)) +
    labs(title = title_text, subtitle = subtitle_text,
         x = "Longitude", y = "Latitude") +
    theme_minimal(base_size = 11) +
    theme(legend.position  = "right",
          panel.grid.major = element_line(color = "grey88", linewidth = 0.3),
          plot.title       = element_text(face = "bold", size = 12),
          plot.subtitle    = element_text(size = 9, color = "grey40"))
}

fig3a <- make_india_map(
  pred_india_df,
  expression(bold("A") ~ " Global model — India"),
  "Crosses = thinned occurrence localities"
)

fig3b <- make_india_map(
  pred_nat_india_df,
  expression(bold("B") ~ " Native range model — India"),
  "Projection of Afrotropical-trained model onto India"
)

fig3 <- grid.arrange(fig3a, fig3b, ncol = 2,
                     top = grid::textGrob(
                       expression(italic("Zaprionus indianus") ~
                                    "— Habitat Suitability in India"),
                       gp = grid::gpar(fontsize = 14, fontface = "bold")))

ggsave("C:/Users/Avirup/OneDrive/Desktop/Zap_SDM/outputs/figures/main/Fig3_india_suitability.pdf",
       fig3, width = 12, height = 7, dpi = 300)
ggsave("C:/Users/Avirup/OneDrive/Desktop/Zap_SDM/outputs/figures/main/Fig3_india_suitability.png",
       fig3, width = 12, height = 7, dpi = 300, bg = "white")

cat("Figure 3 saved.\n")

# =========================================
# MAIN TEXT FIGURE 4
# Binary suitability maps
# Panel A: Global binary | Panel B: India binary
# =========================================

binary_global_df <- as.data.frame(binary_global, xy = TRUE, na.rm = TRUE)
colnames(binary_global_df) <- c("lon", "lat", "suitable")
binary_global_df$suitable <- factor(binary_global_df$suitable,
                                    levels = c(0, 1),
                                    labels = c("Unsuitable", "Suitable"))

binary_india_df <- as.data.frame(binary_india, xy = TRUE, na.rm = TRUE)
colnames(binary_india_df) <- c("lon", "lat", "suitable")
binary_india_df$suitable <- factor(binary_india_df$suitable,
                                   levels = c(0, 1),
                                   labels = c("Unsuitable", "Suitable"))

# Fix 1: Replace geom_raster with geom_tile in both panels
fig4a <- ggplot() +
  geom_tile(data = binary_global_df,
            aes(x = lon, y = lat, fill = suitable)) +
  scale_fill_manual(values = c("Unsuitable" = "#D9D9D9",
                               "Suitable"   = "#2166AC"),
                    name = "Suitability") +
  geom_sf(data = world_sf, fill = NA, color = "white", linewidth = 0.15) +
  geom_sf(data = africa_sf, fill = NA, color = "grey50", linewidth = 0.25) +
  coord_sf(xlim = c(-180, 180), ylim = c(-60, 85)) +
  labs(title    = expression(bold("A") ~ " Global binary suitability"),
       subtitle = paste0("MaxTSS threshold = ", round(threshold_final, 3)),
       x = NULL, y = NULL) +
  theme_minimal(base_size = 11) +
  theme(legend.position  = "bottom",
        panel.grid.major = element_line(color = "grey88", linewidth = 0.3),
        plot.title       = element_text(face = "bold", size = 12),
        plot.subtitle    = element_text(size = 9, color = "grey40"))

fig4b <- ggplot() +
  geom_tile(data = binary_india_df,
            aes(x = lon, y = lat, fill = suitable)) +
  scale_fill_manual(values = c("Unsuitable" = "#D9D9D9",
                               "Suitable"   = "#2166AC"),
                    name = "Suitability") +
  geom_sf(data = india_sf, fill = NA, color = "grey40", linewidth = 0.3) +
  coord_sf(xlim = c(68, 97), ylim = c(8, 37)) +
  labs(title    = expression(bold("B") ~ " Binary suitability — India"),
       subtitle = "State boundaries shown",
       x = "Longitude", y = "Latitude") +
  theme_minimal(base_size = 11) +
  theme(legend.position  = "bottom",
        panel.grid.major = element_line(color = "grey88", linewidth = 0.3),
        plot.title       = element_text(face = "bold", size = 12),
        plot.subtitle    = element_text(size = 9, color = "grey40"))

# Fix 2: grid.arrange with explicit title string, not ...
fig4 <- grid.arrange(
  fig4a, fig4b,
  ncol   = 2,
  widths = c(1.6, 1),
  top    = grid::textGrob(
    expression(italic("Zaprionus indianus") ~
                 "— Binary Habitat Suitability (MaxTSS threshold)"),
    gp = grid::gpar(fontsize = 14, fontface = "bold")
  )
)

ggsave("C:/Users/Avirup/OneDrive/Desktop/Zap_SDM/outputs/figures/main/Fig4_binary_maps.png",
       fig4, width = 14, height = 7, dpi = 300, bg = "white")

# =========================================
# MAIN TEXT FIGURE 5
# Ecospat niche dynamics
# Panel A: Niche overlap plot
# Panel B: Equivalency test
# Panel C: Similarity test
# =========================================

pdf("C:/Users/Avirup/OneDrive/Desktop/Zap_SDM/outputs/figures/main/Fig5_niche_dynamics.pdf",
    width = 14, height = 5)
par(mfrow = c(1, 3), mar = c(4, 4, 3, 2))

# Panel A: Niche dynamics plot
ecospat.plot.niche.dyn(grid.clim.nat, grid.clim.inv,
                       quant      = 0.25,
                       interest   = 2,
                       title      = "",
                       name.axis1 = "PC1",
                       name.axis2 = "PC2")
ecospat.shift.centroids(scores.sp.nat, scores.sp.inv,
                        scores.clim.nat, scores.clim.inv,
                        col = "black")
title(main = expression(bold("A") ~ " Niche dynamics"),
      cex.main = 1.1)
legend("topright",
       legend = c("Native niche", "Invasive niche",
                  "Expansion", "Stability", "Unfilling"),
       fill   = c("#2ecc71", "#e74c3c",
                  "#e74c3c", "#2ecc71", "#3498db"),
       bty    = "n", cex = 0.8)

# Panel B: Equivalency test
ecospat.plot.overlap.test(eq.test, "D",
                          "B  Niche equivalency test")

# Panel C: Similarity test
ecospat.plot.overlap.test(sim.test, "D",
                          "C  Niche similarity test")

dev.off()

# Also save as PNG
png("C:/Users/Avirup/OneDrive/Desktop/Zap_SDM/outputs/figures/main/Fig5_niche_dynamics.png",
    width = 14, height = 5, units = "in", res = 300, bg = "white")
par(mfrow = c(1, 3), mar = c(4, 4, 3, 2))

ecospat.plot.niche.dyn(grid.clim.nat, grid.clim.inv,
                       quant      = 0.25,
                       interest   = 2,
                       title      = "",
                       name.axis1 = "PC1",
                       name.axis2 = "PC2")
# After ecospat.plot.niche.dyn() call, add:
legend("bottomright",
       legend = c("Native niche", "Invasive niche",
                  "Stability", "Expansion", "Unfilling"),
       fill   = c(NA, NA, "blue", "red", "grey"),
       border = c("green", "red", NA, NA, NA),
       bty    = "n", cex = 0.75)
ecospat.shift.centroids(scores.sp.nat, scores.sp.inv,
                        scores.clim.nat, scores.clim.inv,
                        col = "black")
title(main = expression(bold("A") ~ " Niche dynamics"), cex.main = 1.1)

ecospat.plot.overlap.test(eq.test, "D", "B  Niche equivalency test")
ecospat.plot.overlap.test(sim.test, "D", "C  Niche similarity test")

dev.off()
cat("Figure 5 saved.\n")

# =========================================
# MAIN TEXT FIGURE 6 — MESS (revised)
# =========================================

mess.df  <- as.data.frame(mess.object)
colname5 <- colnames(mess.df)[5]

p1 <- ggplot() +
  geom_point(data = mess.df,
             aes(x = lon, y = lat, color = MESS),
             size = 0.8, alpha = 0.7) +
  scale_color_gradient2(low      = "#d73027",
                        mid      = "white",
                        high     = "#4575b4",
                        midpoint = 0,
                        name     = "MESS",
                        limits   = c(-231, 62),
                        oob      = scales::squish) +
  geom_sf(data = world_sf,
          fill = NA, color = "grey40", linewidth = 0.2) +
  geom_sf(data = africa_sf,
          fill = "grey85", color = "grey40",
          linewidth = 0.3, alpha = 0.5) +
  annotate("text", x = 20, y = 10,
           label    = "Native range",
           size     = 3.5,
           color    = "grey20",
           fontface = "bold.italic") +
  coord_sf(xlim = c(-180, 180), ylim = c(-58, 80)) +
  labs(title    = "Multivariate Environmental Similarity Surface (MESS)",
       subtitle = paste0("Blue = similar to native range  |  ",
                         "Red = novel conditions  |  ",
                         "54.6% of invaded range is novel (MESS < 0)"),
       x = NULL, y = NULL) +
  theme_minimal(base_size = 11) +
  theme(legend.position  = "right",
        panel.grid.major = element_line(color = "grey90", linewidth = 0.3),
        plot.title       = element_text(face = "bold", size = 12),
        plot.subtitle    = element_text(size = 9, color = "grey40"))

p2 <- ggplot() +
  geom_point(data = mess.df,
             aes(x = lon, y = lat, color = MESSw),
             size = 0.8, alpha = 0.7) +
  scale_color_gradient2(low      = "#d73027",
                        mid      = "white",
                        high     = "#4575b4",
                        midpoint = 0,
                        name     = "MESSw") +
  geom_sf(data = world_sf,
          fill = NA, color = "grey40", linewidth = 0.2) +
  geom_sf(data = africa_sf,
          fill = "grey85", color = "grey40",
          linewidth = 0.3, alpha = 0.5) +
  annotate("text", x = 20, y = 10,
           label    = "Native range",
           size     = 3.5,
           color    = "grey20",
           fontface = "bold.italic") +
  coord_sf(xlim = c(-180, 180), ylim = c(-58, 80)) +
  labs(title = "Weighted MESS (MESSw)",
       x = NULL, y = NULL) +
  theme_minimal(base_size = 11) +
  theme(legend.position  = "right",
        panel.grid.major = element_line(color = "grey90", linewidth = 0.3),
        plot.title       = element_text(face = "bold", size = 12))

p3 <- ggplot() +
  geom_point(data = mess.df,
             aes(x = lon, y = lat,
                 color = .data[[colname5]]),
             size = 0.8, alpha = 0.7) +
  scale_color_gradient(low    = "#ffffcc",
                       high   = "#d73027",
                       name   = "# Novel\nVariables",
                       limits = c(0, 8),
                       breaks = c(0, 2, 4, 6, 8)) +
  geom_sf(data = world_sf,
          fill = NA, color = "grey40", linewidth = 0.2) +
  geom_sf(data = africa_sf,
          fill = "grey85", color = "grey40",
          linewidth = 0.3, alpha = 0.5) +
  annotate("text", x = 20, y = 10,
           label    = "Native range",
           size     = 3.5,
           color    = "grey20",
           fontface = "bold.italic") +
  coord_sf(xlim = c(-180, 180), ylim = c(-58, 80)) +
  labs(title    = "Number of variables with novel conditions",
       subtitle = "Higher values indicate greater extrapolation risk",
       x = NULL, y = NULL) +
  theme_minimal(base_size = 11) +
  theme(legend.position  = "right",
        panel.grid.major = element_line(color = "grey90", linewidth = 0.3),
        plot.title       = element_text(face = "bold", size = 12),
        plot.subtitle    = element_text(size = 9, color = "grey40"))

# Combine and save
mess_combined <- grid.arrange(
  p1, p2, p3,
  ncol = 1,
  top  = grid::textGrob(
    expression(italic("Zaprionus indianus") ~
                 "— Climatic Novelty in the Invaded Range"),
    gp = grid::gpar(fontsize = 14, fontface = "bold")
  )
)

ggsave("C:/Users/Avirup/OneDrive/Desktop/Zap_SDM/outputs/figures/main/Fig6_MESS.png",
       plot   = mess_combined,
       width  = 12,
       height = 14,
       dpi    = 300,
       bg     = "white")

ggsave("C:/Users/Avirup/OneDrive/Desktop/Zap_SDM/outputs/figures/main/Fig6_MESS.pdf",
       plot   = mess_combined,
       width  = 12,
       height = 14,
       dpi    = 300)

cat("Figure 6 saved.\n")

# =========================================
# SUPPLEMENTARY FIGURE S1
# Pearson correlation matrix — all 19 variables
# =========================================

env_vals_all <- extract(env,
                        as.matrix(occ_thin[, c("lon", "lat")]))[, names(env)]
colnames(env_vals_all) <- names(env)
cor_mat_all <- cor(env_vals_all, use = "complete.obs")

pdf("C:/Users/Avirup/OneDrive/Desktop/Zap_SDM/outputs/figures/supplementary/FigS1_correlation_matrix.pdf",
    width = 10, height = 9)
corrplot(cor_mat_all,
         method   = "color",
         type     = "upper",
         tl.cex   = 0.85,
         tl.col   = "black",
         addCoef.col = "black",
         number.cex  = 0.55,
         col      = colorRampPalette(c("#d73027","white","#4575b4"))(200),
         title    = "Pearson correlation — 19 bioclimatic variables",
         mar      = c(0, 0, 2, 0))
dev.off()

png("C:/Users/Avirup/OneDrive/Desktop/Zap_SDM/outputs/figures/supplementary/FigS1_correlation_matrix.png",
    width = 10, height = 9, units = "in", res = 300, bg = "white")
corrplot(cor_mat_all,
         method      = "color",
         type        = "upper",
         tl.cex      = 0.85,
         tl.col      = "black",
         addCoef.col = "black",
         number.cex  = 0.55,
         col         = colorRampPalette(c("#d73027","white","#4575b4"))(200),
         title       = "Pearson correlation — 19 bioclimatic variables",
         mar         = c(0, 0, 2, 0))
dev.off()
cat("Figure S1 saved.\n")

# =========================================
# SUPPLEMENTARY FIGURE S2
# VIF values — retained variables
# =========================================

vif_df <- data.frame(
  variable = final_vars,
  VIF      = c(3.869736, 3.139302, 3.903550, 3.109993,
               1.889089, 3.188780, 2.405857, 3.481678)
)

figS2 <- ggplot(vif_df, aes(x = reorder(variable, VIF), y = VIF)) +
  geom_col(fill = "#4393c3", color = "grey30", width = 0.6) +
  geom_hline(yintercept = 10, linetype = "dashed",
             color = "red", linewidth = 0.8) +
  geom_hline(yintercept = 5,  linetype = "dotted",
             color = "orange", linewidth = 0.7) +
  geom_text(aes(label = round(VIF, 2)), hjust = -0.2,
            size = 3.5, color = "grey20") +
  coord_flip(ylim = c(0, 12)) +
  annotate("text", x = 1, y = 10.2, label = "VIF = 10 threshold",
           color = "red", size = 3, hjust = 0) +
  labs(title    = "Variance Inflation Factor — retained variables",
       subtitle = "All variables VIF < 10; dashed line = exclusion threshold",
       x = "Bioclimatic variable", y = "VIF") +
  theme_minimal(base_size = 11) +
  theme(plot.title    = element_text(face = "bold", size = 12),
        plot.subtitle = element_text(size = 9, color = "grey40"))

ggsave("C:/Users/Avirup/OneDrive/Desktop/Zap_SDM/outputs/figures/supplementary/FigS2_VIF.pdf",
       figS2, width = 7, height = 5, dpi = 300)
ggsave("C:/Users/Avirup/OneDrive/Desktop/Zap_SDM/outputs/figures/supplementary/FigS2_VIF.png",
       figS2, width = 7, height = 5, dpi = 300, bg = "white")
cat("Figure S2 saved.\n")

# =========================================
# SUPPLEMENTARY FIGURE S3
# Response curves — global vs native model
# Temperature variables divided by 10 for correct axis display
# =========================================

# Variables to divide by 10 for display
temp_vars <- c("bio5", "bio9", "bio10", "bio11")

var_labels <- c(
  bio5  = "Max Temp Warmest Month (°C)",
  bio9  = "Mean Temp Driest Quarter (°C)",
  bio10 = "Mean Temp Warmest Quarter (°C)",
  bio11 = "Mean Temp Coldest Quarter (°C)",
  bio12 = "Annual Precipitation (mm)",
  bio13 = "Precip Wettest Month (mm)",
  bio18 = "Precip Warmest Quarter (mm)",
  bio19 = "Precip Coldest Quarter (mm)"
)

plot_response <- function() {
  par(mfrow = c(3, 3), mar = c(4, 4, 3, 1))
  
  for (v in final_vars) {
    
    x <- seq(min(dat_full[[v]], na.rm = TRUE),
             max(dat_full[[v]], na.rm = TRUE),
             length.out = 200)
    
    newdata           <- as.data.frame(
      matrix(rep(colMeans(dat_full[, final_vars], na.rm = TRUE),
                 each = 200), nrow = 200))
    colnames(newdata) <- final_vars
    newdata[[v]]      <- x
    
    y_full <- predict(model_all, newdata, type = "cloglog")
    y_nat  <- predict(model_nat, newdata, type = "cloglog")
    
    # Divide temperature x-axis by 10 for correct display
    x_display <- if (v %in% temp_vars) x / 10 else x
    
    # Divide rug values by 10 for temperature variables
    rug_vals <- dat_full[dat_full$pa == 1, v]
    rug_display <- if (v %in% temp_vars) rug_vals / 10 else rug_vals
    
    plot(x_display, y_full,
         type     = "l", lwd = 2.5, col = "#e74c3c",
         ylim     = c(0, max(y_full, y_nat, na.rm = TRUE) * 1.1),
         xlab     = var_labels[v],
         ylab     = "Habitat suitability",
         main     = v,
         cex.main = 1.1, cex.lab = 0.9,
         bty      = "l")
    lines(x_display, y_nat,
          col = "#2166ac", lwd = 2.5, lty = 2)
    rug(rug_display,
        col      = "grey40",
        ticksize = 0.03)
  }
  
  # Legend panel
  plot.new()
  legend("center",
         legend = c("Global model (all occurrences)",
                    "Native range model (Africa only)"),
         col    = c("#e74c3c", "#2166ac"),
         lwd    = 2.5,
         lty    = c(1, 2),
         bty    = "n",
         cex    = 1.0,
         title  = expression(bold("Model")))
}

# Save PDF
pdf("C:/Users/Avirup/OneDrive/Desktop/Zap_SDM/outputs/figures/supplementary/FigS3_response_curves.pdf",
    width = 12, height = 10)
plot_response()
dev.off()

# Save PNG
png("C:/Users/Avirup/OneDrive/Desktop/Zap_SDM/outputs/figures/supplementary/FigS3_response_curves.png",
    width = 12, height = 10, units = "in", res = 300, bg = "white")
plot_response()
dev.off()

cat("Figure S3 saved.\n")

# =========================================
# SUPPLEMENTARY FIGURE S4
# PCA contribution plot (ecospat)
# =========================================

# Pre-compute title string once
pca_title <- paste0(
  "PCA — bioclimatic variable contributions\n",
  "PC1: ",        round(pca.env$eig[1] / sum(pca.env$eig) * 100, 1),
  "%  |  PC2: ",  round(pca.env$eig[2] / sum(pca.env$eig) * 100, 1),
  "%  |  Cumulative: ",
  round(sum(pca.env$eig[1:2]) / sum(pca.env$eig) * 100, 1), "%"
)

plot_pca <- function() {
  # Increase top margin to make room for two-line custom title
  par(mar = c(5, 4, 5, 2))
  
  # Draw ecospat plot
  ecospat.plot.contrib(contrib = pca.env$co, eigen = pca.env$eig)
  
  # Blank out ecospat's default "correlation circle" title
  title(main = "", line = 3)
  
  # Add custom title on top
  title(main = pca_title, cex.main = 1.0, line = 3)
}

# Save PDF
pdf("C:/Users/Avirup/OneDrive/Desktop/Zap_SDM/outputs/figures/supplementary/FigS4_PCA_contribution.pdf",
    width = 7, height = 7.5)
plot_pca()
dev.off()

# Save PNG
png("C:/Users/Avirup/OneDrive/Desktop/Zap_SDM/outputs/figures/supplementary/FigS4_PCA_contribution.png",
    width = 7, height = 7.5, units = "in", res = 300, bg = "white")
plot_pca()
dev.off()

cat("Figure S4 saved.\n")

# =========================================
# SUPPLEMENTARY FIGURE S5
# Univariate niche plots
# =========================================

nat_pres_plot <- nat[nat$species == 1, ]
inv_pres_plot <- inv[inv$species == 1, ]

pdf("C:/Users/Avirup/OneDrive/Desktop/Zap_SDM/outputs/figures/supplementary/FigS5_univariate_niche.pdf",
    width = 12, height = 10)
par(mfrow = c(3, 3), mar = c(4, 4, 3, 1))

for (v in env_cols) {
  d_nat_bg <- density(nat[[v]],          na.rm = TRUE)
  d_inv_bg <- density(inv[[v]],          na.rm = TRUE)
  d_nat    <- density(nat_pres_plot[[v]], na.rm = TRUE)
  d_inv    <- density(inv_pres_plot[[v]], na.rm = TRUE)
  
  ylim_max <- max(d_nat_bg$y, d_inv_bg$y,
                  d_nat$y,    d_inv$y, na.rm = TRUE) * 1.1
  
  plot(d_nat_bg,
       col  = "grey50", lwd  = 1.5, lty  = 1,
       main = var_labels[v],
       xlab = var_labels[v], ylab = "Density",
       ylim = c(0, ylim_max),
       cex.main = 0.9, bty = "l")
  lines(d_inv_bg, col = "grey50", lwd = 1.5, lty = 2)
  lines(d_nat,    col = col_native,   lwd = 2.5)
  lines(d_inv,    col = col_invasive, lwd = 2.5)
  
  x_common <- seq(min(d_nat$x, d_inv$x),
                  max(d_nat$x, d_inv$x), length = 500)
  y_nat_i  <- approx(d_nat$x, d_nat$y, xout = x_common)$y
  y_inv_i  <- approx(d_inv$x, d_inv$y, xout = x_common)$y
  
  polygon(c(x_common, rev(x_common)),
          c(pmin(y_nat_i, y_inv_i, na.rm = TRUE),
            rep(0, length(x_common))),
          col    = adjustcolor("purple", 0.35),
          border = NA)
}

plot.new()
legend("center",
       legend = c("Native climate space", "Invasive climate space",
                  "Native niche", "Invasive niche", "Overlap"),
       col    = c("grey50", "grey50",
                  col_native, col_invasive, "purple"),
       lwd    = c(1.5, 1.5, 2.5, 2.5, 6),
       lty    = c(1, 2, 1, 1, 1),
       bty    = "n", cex = 0.95,
       title  = expression(bold("Legend")))

dev.off()

png("C:/Users/Avirup/OneDrive/Desktop/Zap_SDM/outputs/figures/supplementary/FigS5_univariate_niche.png",
    width = 12, height = 10, units = "in", res = 300, bg = "white")
par(mfrow = c(3, 3), mar = c(4, 4, 3, 1))
for (v in env_cols) {
  d_nat_bg <- density(nat[[v]], na.rm = TRUE)
  d_inv_bg <- density(inv[[v]], na.rm = TRUE)
  d_nat    <- density(nat_pres_plot[[v]], na.rm = TRUE)
  d_inv    <- density(inv_pres_plot[[v]], na.rm = TRUE)
  ylim_max <- max(d_nat_bg$y, d_inv_bg$y,
                  d_nat$y, d_inv$y, na.rm = TRUE) * 1.1
  plot(d_nat_bg, col = "grey50", lwd = 1.5, lty = 1,
       main = var_labels[v], xlab = var_labels[v],
       ylab = "Density", ylim = c(0, ylim_max),
       cex.main = 0.9, bty = "l")
  lines(d_inv_bg, col = "grey50", lwd = 1.5, lty = 2)
  lines(d_nat, col = col_native,   lwd = 2.5)
  lines(d_inv, col = col_invasive, lwd = 2.5)
  x_common <- seq(min(d_nat$x, d_inv$x),
                  max(d_nat$x, d_inv$x), length = 500)
  y_nat_i <- approx(d_nat$x, d_nat$y, xout = x_common)$y
  y_inv_i <- approx(d_inv$x, d_inv$y, xout = x_common)$y
  polygon(c(x_common, rev(x_common)),
          c(pmin(y_nat_i, y_inv_i, na.rm = TRUE),
            rep(0, length(x_common))),
          col = adjustcolor("purple", 0.35), border = NA)
}
plot.new()
legend("center",
       legend = c("Native climate space", "Invasive climate space",
                  "Native niche", "Invasive niche", "Overlap"),
       col  = c("grey50", "grey50", col_native, col_invasive, "purple"),
       lwd  = c(1.5, 1.5, 2.5, 2.5, 6),
       lty  = c(1, 2, 1, 1, 1),
       bty  = "n", cex = 0.95)
dev.off()
cat("Figure S5 saved.\n")

# =========================================
# SUPPLEMENTARY FIGURE S6
# Model evaluation — ROC curves
# =========================================

pred_vals_full <- as.numeric(predict(model_all, dat_full[, final_vars],
                                     type = "cloglog"))
roc_full_plot  <- pROC::roc(dat_full$pa, pred_vals_full, quiet = TRUE)
roc_test_plot  <- pROC::roc(test_data$pa, pred_test, quiet = TRUE)

roc_full_df <- data.frame(
  FPR   = 1 - roc_full_plot$specificities,
  TPR   = roc_full_plot$sensitivities,
  model = paste0("Training (AUC = ",
                 round(pROC::auc(roc_full_plot), 3), ")")
)

roc_test_df <- data.frame(
  FPR   = 1 - roc_test_plot$specificities,
  TPR   = roc_test_plot$sensitivities,
  model = paste0("Test 70/30 (AUC = ",
                 round(pROC::auc(roc_test_plot), 3), ")")
)

roc_df <- rbind(roc_full_df, roc_test_df)

figS6 <- ggplot(roc_df, aes(x = FPR, y = TPR, color = model)) +
  geom_line(linewidth = 1.2) +
  geom_abline(slope = 1, intercept = 0,
              linetype = "dashed", color = "grey60") +
  scale_color_manual(values = c("#e74c3c", "#2166ac"),
                     name   = "Model") +
  annotate("text", x = 0.65, y = 0.15,
           label = paste0("TSS (training) = ", round(TSS_full, 3),
                          "\nTSS (test) = ",   round(TSS_test, 3)),
           size = 3.5, color = "grey20", hjust = 0) +
  labs(title    = "ROC curves — model evaluation",
       subtitle = expression(italic("Zaprionus indianus") ~
                               "global distribution model"),
       x = "False positive rate (1 - Specificity)",
       y = "True positive rate (Sensitivity)") +
  theme_minimal(base_size = 11) +
  theme(legend.position  = "bottom",
        panel.grid.major = element_line(color = "grey88", linewidth = 0.3),
        plot.title       = element_text(face = "bold", size = 12),
        plot.subtitle    = element_text(size = 9, color = "grey40"))

ggsave("C:/Users/Avirup/OneDrive/Desktop/Zap_SDM/outputs/figures/supplementary/FigS6_ROC_curves.pdf",
       figS6, width = 7, height = 6, dpi = 300)
ggsave("C:/Users/Avirup/OneDrive/Desktop/Zap_SDM/outputs/figures/supplementary/FigS6_ROC_curves.png",
       figS6, width = 7, height = 6, dpi = 300, bg = "white")
cat("Figure S6 saved.\n")

# =========================================
# SUMMARY OF ALL FIGURES
# =========================================

cat("\n========================================\n")
cat("ALL FIGURES SAVED\n")
cat("========================================\n")
cat("\nMAIN TEXT:\n")
cat("Fig 1 — Global occurrence map\n")
cat("Fig 2 — Global habitat suitability (continuous)\n")
cat("Fig 3 — India habitat suitability\n")
cat("Fig 4 — Binary suitability maps\n")
cat("Fig 5 — Niche dynamics (ecospat)\n")
cat("Fig 6 — MESS analysis\n")
cat("\nSUPPLEMENTARY:\n")
cat("Fig S1 — Pearson correlation matrix (19 variables)\n")
cat("Fig S2 — VIF values (retained variables)\n")
cat("Fig S3 — Response curves (global vs native model)\n")
cat("Fig S4 — PCA contribution plot\n")
cat("Fig S5 — Univariate niche plots\n")
cat("Fig S6 — ROC curves and evaluation metrics\n")
cat("\nLocation: C:/Users/Avirup/OneDrive/Desktop/Zap_SDM/outputs/figures/\n")