R-4.5.2 
SOFTWARE REQUIREMENTS
These scripts were developed and tested in R version 4.2.0+. 
The following libraries are required:
- terra (Spatial data handling)
- biomod2 (Species Distribution Modeling)
- humboldt (Niche overlap and divergence tests)
- dplyr, tidyr (Data manipulation)

# ==============================================================================
# HUMBOLDT ANALYTICAL PROTOCOL: NICHE EVOLUTION & DIVERGENCE (FINAL VERSION)
# ==============================================================================
# Purpose: Quantifying niche overlap and divergence between L. canadensis (Lc) 
# and L. rufus (Lr) in E-space using PCA-based transformations.
# Rationale: This protocol avoids RAM exhaustion by using tabular environmental 
# background data and applies 300 DPI resolution for publication-quality outputs.
# ==============================================================================

library(humboldt)
library(dplyr)

# --- 1. DIRECTORY SETUP ---
# Defining paths for DATA structure
occ_dir <- "C:/DATA/Humboldt_data"
out_dir <- "C:/OUTPUTS"

if (!dir.exists(out_dir)) dir.create(out_dir, recursive = TRUE)

# ------------------------------------------------------------------------------
# STEP 1 & 2: LOADING PRE-PROCESSED CSV DATA
# ------------------------------------------------------------------------------
message(">>> Loading high-precision occurrence and environmental background data...")

sp1  <- read.csv(file.path(occ_dir, "Humboldt_sp1_Lc.csv"))
sp2  <- read.csv(file.path(occ_dir, "Humboldt_sp2_Lr.csv"))
env1 <- read.csv(file.path(occ_dir, "Humboldt_env1_Lc.csv"))
env2 <- read.csv(file.path(occ_dir, "Humboldt_env2_Lr.csv"))

# Define column range for environmental variables (starting from column 3)
e_vars_cols <- 3:ncol(env1)

# ------------------------------------------------------------------------------
# STEP 3: G2E TRANSFORMATION (GEOGRAPHIC TO ENVIRONMENTAL SPACE)
# ------------------------------------------------------------------------------
# Rationale: Transforming G-space to E-space via PCA to account for environmental 
# availability and reduce dimensionality.
message(">>> Executing G2E transformation (PCA-based E-space creation)...")

zz <- humboldt.g2e(env1 = env1, env2 = env2, sp1 = sp1, sp2 = sp2, 
                   reduce.env = 2, reductype = "PCA", 
                   non.analogous.environments = "NO", 
                   env.trim = TRUE, e.var = e_vars_cols, col.env = e_vars_cols, 
                   trim.buffer.sp1 = 500, trim.buffer.sp2 = 500, 
                   rarefy.dist = 50, rarefy.units = "km", env.reso = 0.041666667, 
                   kern.smooth = 1, R = 100, run.silent = FALSE)

# Extract scores for the first two principal components
scores_env12 <- rbind(zz$scores.env1[1:2], zz$scores.env2[1:2])
scores_env1  <- zz$scores.env1[1:2]
scores_env2  <- zz$scores.env2[1:2]
scores_sp1   <- zz$scores.sp1[1:2]
scores_sp2   <- zz$scores.sp2[1:2]

# ------------------------------------------------------------------------------
# STEP 4 & 5: ENVIRONMENTAL GRIDDING AND PRIMARY STATISTICS (NOT & PNTI)
# ------------------------------------------------------------------------------
message(">>> Creating environmental grids and calculating niche statistics...")

z_sp1 <- humboldt.grid.espace(scores_env12, scores_env1, scores_sp1, kern.smooth=1, R=100)
z_sp2 <- humboldt.grid.espace(scores_env12, scores_env2, scores_sp2, kern.smooth=1, R=100)

# Calculate niche similarity (Schoener's D) and Equivalency
sim_correct <- humboldt.niche.similarity(z_sp1, z_sp2, correct.env=TRUE)
eq_test     <- humboldt.equivalence.stat(z_sp1, z_sp2, rep=100, kern.smooth=1, ncores=2)

# Calculate Potential Niche Truncation Index (PNTI)
# High PNTI indicates niche limits are constrained by available environmental space
pnt_lc <- humboldt.pnt.index(scores_env12, scores_env1, scores_sp1, kern.smooth=1, R=100)
pnt_lr <- humboldt.pnt.index(scores_env12, scores_env2, scores_sp2, kern.smooth=1, R=100)

# ------------------------------------------------------------------------------
# STEP 6: NICHE DIVERGENCE TEST (NDT / BACKGROUND TEST)
# ------------------------------------------------------------------------------
# Rationale: Testing if observed niche differences are significant relative to 
# the differences in available environmental backgrounds.
message(">>> Running NDT (Background) tests (100 iterations per direction)...")

bg_1to2 <- humboldt.background.stat(g2e = zz, rep = 100, sim.dir = 1, 
                                    env.reso = 0.041666667, kern.smooth = 1, 
                                    correct.env = TRUE, R = 100, run.silent.bak = FALSE)

bg_2to1 <- humboldt.background.stat(g2e = zz, rep = 100, sim.dir = 2, 
                                    env.reso = 0.041666667, kern.smooth = 1, 
                                    correct.env = TRUE, R = 100, run.silent.bak = FALSE)

# ------------------------------------------------------------------------------
# STEP 7-9: VISUALIZATION & EXPORT (HIGH-RESOLUTION PDF & PNG)
# ------------------------------------------------------------------------------
message(">>> Generating publication-quality graphics (300 DPI)...")

z_env1 <- humboldt.grid.espace(scores_env12, scores_env1, scores_env1, kern.smooth=1, R=100)
z_env2 <- humboldt.grid.espace(scores_env12, scores_env2, scores_env2, kern.smooth=1, R=100)
ee     <- humboldt.espace.correction(Z.env1 = z_env1, Z.env2 = z_env2, Z.sp1 = z_sp1, Z.sp2 = z_sp2)

# Prepare custom difference plot (Lc=Blue, Lr=Red)
ee_custom <- ee
temp_s1 <- ee_custom$s1
ee_custom$s1 <- ee_custom$s2
ee_custom$s2 <- temp_s1

# Function to save plots as both PDF and PNG
export_niche_plot <- function(name, plot_func, width=8, height=8) {
  # PDF
  pdf(file.path(out_dir, paste0(name, ".pdf")), width=width, height=height)
  plot_func()
  dev.off()
  # PNG (300 DPI)
  png(file.path(out_dir, paste0(name, ".png")), width=width*300, height=height*300, res=300)
  plot_func()
  dev.off()
}

export_niche_plot("Humboldt_1_Niche_Lc", function() humboldt.plot.niche(z_sp1, "L. canadensis (Lc)", "PC1", "PC2"))
export_niche_plot("Humboldt_2_Niche_Lr", function() humboldt.plot.niche(z_sp2, "L. rufus (Lr)", "PC1", "PC2"))
export_niche_plot("Humboldt_3_Espace_Difference", function() humboldt.plot.espace.diff(ee_custom, correct.env = TRUE, type = "species"))

# Export Density Plots for NDT
png(file.path(out_dir, "Humboldt_4_NDT_Density.png"), width=3000, height=1500, res=300)
par(mfrow=c(1,2))
humboldt.plot.density(bg_1to2, "D", "Background 1 -> 2 (Lc to Lr)")
humboldt.plot.density(bg_2to1, "D", "Background 2 -> 1 (Lr to Lc)")
dev.off()

# ------------------------------------------------------------------------------
# FINAL STATISTICAL SUMMARY
# ------------------------------------------------------------------------------
cat("\n=================================================================\n")
cat("                HUMBOLDT ANALYSIS COMPLETE SUMMARY               \n")
cat("=================================================================\n")
cat("1. Schoener's D (Overlap)          :", sim_correct$D, "\n")
cat("2. Equivalency (NOT) P-Value       :", eq_test$p.D, "\n")
cat("3. NDT 1->2 (Lc->Lr) P-Value       :", bg_1to2$p.D, "\n")
cat("4. NDT 2->1 (Lr->Lc) P-Value       :", bg_2to1$p.D, "\n")
cat("5. Lc PNTI (Truncation Index)      :", pnt_lc$pnt.index, "\n")
cat("6. Lr PNTI (Truncation Index)      :", pnt_lr$pnt.index, "\n")
cat("-----------------------------------------------------------------\n")
message(">>> All PDF/PNG outputs successfully exported to: ", out_dir)







# ------------------------------------------------------------------------------
# STEP 10: EXPORT FINAL STATISTICAL REPORT TO TEXT FILE
# ------------------------------------------------------------------------------
message(">>> Exporting final statistical summary to text file...")

# Define output file path
report_path <- file.path(out_dir, "Humboldt_Final_Statistical_Report.txt")

# Start sinking output to the file
sink(report_path)

cat("=================================================================\n")
cat("          HUMBOLDT NICHE ANALYSIS: STATISTICAL SUMMARY           \n")
cat("                Species: L. canadensis vs L. rufus               \n")
cat("=================================================================\n\n")

cat("1. NICHE OVERLAP (Schoener's D)\n")
cat("-----------------------------------------------------------------\n")
cat("Observed Schoener's D: ", round(sim_correct$D, 4), "\n")
cat("Interpretation: 0 = No overlap, 1 = Complete overlap\n\n")

cat("2. NICHE EQUIVALENCY TEST (NOT)\n")
cat("-----------------------------------------------------------------\n")
cat("Null Hypothesis: The niches are identical.\n")
cat("P-value (D-metric)   : ", round(eq_test$p.D, 4), "\n")
cat("Decision             : ", ifelse(eq_test$p.D < 0.05, 
                                      "REJECT Null (Niches are significantly different)", 
                                      "FAIL TO REJECT Null (Niches are equivalent)"), "\n\n")

cat("3. NICHE DIVERGENCE TEST (NDT / BACKGROUND TEST)\n")
cat("-----------------------------------------------------------------\n")
cat("Null Hypothesis: Niche differences are caused by environmental availability.\n")
cat("Direction 1 (Lc -> Lr) P-value: ", round(bg_1to2$p.D, 4), "\n")
cat("Direction 2 (Lr -> Lc) P-value: ", round(bg_2to1$p.D, 4), "\n")
cat("Decision (Lc -> Lr)  : ", ifelse(bg_1to2$p.D < 0.05, "Niche Divergence", "Niche Conservatism"), "\n")
cat("Decision (Lr -> Lc)  : ", ifelse(bg_2to1$p.D < 0.05, "Niche Divergence", "Niche Conservatism"), "\n\n")

cat("4. POTENTIAL NICHE TRUNCATION INDEX (PNTI)\n")
cat("-----------------------------------------------------------------\n")
cat("Rationale: Measures if the niche is cut off by available environmental limits.\n")
cat("PNTI for L. canadensis (Lc): ", round(pnt_lc$pnt.index, 4), "\n")
cat("PNTI for L. rufus (Lr)     : ", round(pnt_lr$pnt.index, 4), "\n")
cat("Note: Values near 1.0 indicate strong truncation (abiotic constraints).\n\n")

cat("=================================================================\n")
cat("Report Generated on: ", format(Sys.time(), "%Y-%m-%d %H:%M:%S"), "\n")
cat("Protocol: Humboldt Manual G2E (Sari et al., 2026)\n")
cat("=================================================================\n")

# Stop sinking
sink()

message(">>> Final report muzzled and saved to: ", report_path)


