1). Descarga y preparación de datos de presencia

Primero se usó el paquete rgbif para descargar registros georreferenciados del hongo Fusarium oxysporum f. sp. cubense desde la Global Biodiversity Information Facility (GBIF).
Estos registros representan los puntos de presencia confirmada del patógeno en distintas regiones del mundo.

Luego, mediante CoordinateCleaner, se eliminaron registros erróneos o dudosos (coordenadas en instituciones, ceros, duplicados o centroides de países).
El resultado fueron coordenadas ecológicamente válidas, representadas como objetos espaciales con sistema de referencia WGS84 (EPSG:4326).

Conceptualmente, esta etapa representa el espacio observado de ocurrencias reales, el cual será usado para estimar la probabilidad de presencia del organismo.

# --- INSTALAR Y CARGAR LIBRERÍAS NECESARIAS ---
if (!requireNamespace("rgbif")) install.packages("rgbif")
## Loading required namespace: rgbif
if (!requireNamespace("CoordinateCleaner")) install.packages("CoordinateCleaner")
## Loading required namespace: CoordinateCleaner
if (!requireNamespace("geodata")) install.packages("geodata")
## Loading required namespace: geodata
if (!requireNamespace("raster")) install.packages("raster")
## Loading required namespace: raster
if (!requireNamespace("dismo")) install.packages("dismo")
## Loading required namespace: dismo
if (!requireNamespace("maxnet")) install.packages("maxnet")
## Loading required namespace: maxnet
if (!requireNamespace("ENMeval")) install.packages("ENMeval")
## Loading required namespace: ENMeval
if (!requireNamespace("sf")) install.packages("sf")
## Loading required namespace: sf
if (!requireNamespace("tidyverse")) install.packages("tidyverse")
## Loading required namespace: tidyverse
library(rgbif)
## Warning: package 'rgbif' was built under R version 4.4.3
library(CoordinateCleaner)
## Warning: package 'CoordinateCleaner' was built under R version 4.4.3
library(geodata)
## Warning: package 'geodata' was built under R version 4.4.3
## Loading required package: terra
## Warning: package 'terra' was built under R version 4.4.3
## terra 1.8.80
library(raster)
## Loading required package: sp
## Warning: package 'sp' was built under R version 4.4.3
library(dismo)
## Warning: package 'dismo' was built under R version 4.4.3
library(maxnet)
## Warning: package 'maxnet' was built under R version 4.4.3
library(ENMeval)
## Warning: package 'ENMeval' was built under R version 4.4.3
## This is ENMeval version 2.0.5.2. 
## For worked examples, please consult the vignette: <https://jamiemkass.github.io/ENMeval/articles/ENMeval-2.0-vignette.html>.
library(sf)
## Warning: package 'sf' was built under R version 4.4.3
## Linking to GEOS 3.13.0, GDAL 3.10.1, PROJ 9.5.1; sf_use_s2() is TRUE
library(tidyverse)
## Warning: package 'ggplot2' was built under R version 4.4.3
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   4.0.0     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.1
## ✔ purrr     1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ tidyr::extract() masks raster::extract(), terra::extract()
## ✖ dplyr::filter()  masks stats::filter()
## ✖ dplyr::lag()     masks stats::lag()
## ✖ dplyr::select()  masks raster::select()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
get_occurrences <- function(species_name, limit = 5000) {
  data <- occ_search(
    scientificName = species_name,
    hasCoordinate = TRUE,
    limit = limit
  )
  df <- data$data %>%
    select(decimalLongitude, decimalLatitude, country, year) %>%
    filter(!is.na(decimalLongitude), !is.na(decimalLatitude))
  df$species <- species_name
  df_clean <- clean_coordinates(
    x = df,
    lon = "decimalLongitude",
    lat = "decimalLatitude",
    species = "species",
    tests = c("centroids", "equal", "gbif", "institutions", "zeros")
  ) %>% filter(.summary == TRUE)
  return(df_clean)
}

# Fusarium patógeno
occ_foc <- get_occurrences("Fusarium oxysporum f.sp. cubense")
## Testing coordinate validity
## Flagged 0 records.
## Testing equal lat/lon
## Flagged 0 records.
## Testing zero coordinates
## Flagged 0 records.
## Testing country centroids
## Flagged 0 records.
## Testing GBIF headquarters, flagging records around Copenhagen
## Flagged 0 records.
## Testing biodiversity institutions
## Flagged 1 records.
## Flagged 1 of 5000 records, EQ = 0.
# Hospederos: Musa acuminata (banano) y Musa balbisiana (plátano)
occ_banana <- get_occurrences("Musa acuminata")
## Testing coordinate validity
## Flagged 0 records.
## Testing equal lat/lon
## Flagged 605 records.
## Testing zero coordinates
## Flagged 605 records.
## Testing country centroids
## Flagged 627 records.
## Testing GBIF headquarters, flagging records around Copenhagen
## Flagged 0 records.
## Testing biodiversity institutions
## Flagged 13 records.
## Flagged 640 of 3855 records, EQ = 0.17.
occ_plantain <- get_occurrences("Musa balbisiana")
## Testing coordinate validity
## Flagged 0 records.
## Testing equal lat/lon
## Flagged 45 records.
## Testing zero coordinates
## Flagged 45 records.
## Testing country centroids
## Flagged 52 records.
## Testing GBIF headquarters, flagging records around Copenhagen
## Flagged 0 records.
## Testing biodiversity institutions
## Flagged 10 records.
## Flagged 62 of 584 records, EQ = 0.11.

2). Descarga de variables ambientales (WorldClim v2.1)

Se descargaron las 19 variables bioclimáticas (Bio1–Bio19) desde WorldClim v2.1, que sintetizan información de temperatura y precipitación promedio global.
Estas variables representan dimensiones climáticas relevantes para la ecología de los organismos, por ejemplo:

Se recortaron las capas a la extensión geográfica de interés (–40° a 40° de latitud), generando un RasterStack de predictores ambientales con resolución de 10 minutos (~18 km).

En términos técnicos, este conjunto de capas define el espacio ambiental multivariado (E) dentro del cual el modelo busca estimar la distribución potencial del patógeno.

# --- 2️⃣ DESCARGA DE VARIABLES BIOCLIMÁTICAS DE WORLDCLIM GLOBAL ---
# Usamos las 19 variables BIOCLIM a nivel global (resolución 10')
bio_global <- worldclim_global(var = "bio", res = 10, path = "data")
bio_stack <- stack(bio_global)
names(bio_stack)
##  [1] "wc2.1_10m_bio_1"  "wc2.1_10m_bio_2"  "wc2.1_10m_bio_3"  "wc2.1_10m_bio_4" 
##  [5] "wc2.1_10m_bio_5"  "wc2.1_10m_bio_6"  "wc2.1_10m_bio_7"  "wc2.1_10m_bio_8" 
##  [9] "wc2.1_10m_bio_9"  "wc2.1_10m_bio_10" "wc2.1_10m_bio_11" "wc2.1_10m_bio_12"
## [13] "wc2.1_10m_bio_13" "wc2.1_10m_bio_14" "wc2.1_10m_bio_15" "wc2.1_10m_bio_16"
## [17] "wc2.1_10m_bio_17" "wc2.1_10m_bio_18" "wc2.1_10m_bio_19"
# --- 3️⃣ MAPA EXPLORATORIO DE VARIABLES CLIMÁTICAS ---
plot(bio_stack[[c(1,12,16,17,5,6,7,8,18)]],
     main = names(bio_stack[[c(1,12,16,17,5,6,7,8,18)]]))

3). Generación de puntos de fondo (pseudoausencias)

Dado que los modelos de presencia-ausencia necesitan conocer también dónde no está presente la especie, pero GBIF solo provee datos de presencia, se generaron puntos de fondo (background) aleatorios en el área de estudio mediante randomPoints() de dismo.

Estos puntos no representan ausencias reales, sino una muestra aleatoria del espacio ambiental global, utilizada para contrastar las condiciones climáticas de presencia vs. no-presencia observada.

Matemáticamente, estos puntos aproximan el espacio disponible (M) dentro del marco teórico del nicho de Hutchinson, contra el cual MaxEnt evalúa la probabilidad relativa de ocurrencia.

# --- 4️⃣ CREAR PUNTOS DE FONDO (PSEUDO-AUSENCIAS) ---
set.seed(42)
bg_global <- randomPoints(bio_stack, n = 20000)

4). Modelado de distribución con MaxEnt

El modelo se construyó con la función maxent() del paquete dismo, que implementa el algoritmo de Máxima Entropía (Maximum Entropy Modeling).
Este enfoque estima una función de probabilidad de presencia P(x) sobre el espacio ambiental X que cumple dos condiciones:

  1. Satisface las restricciones empíricas derivadas de los promedios ambientales observados en los puntos de presencia.
    Formalmente:

EP[fi(x)]=μi

  1. donde fi(x)son las variables ambientales y μi sus medias empíricas.

    1. Maximiza la entropía de Shannon del sistema:
    1. H(P)=−∑xP(x)log⁡P(x)
    es decir, el modelo más uniforme posible (menos sesgado) sujeto a las restricciones observadas.

El resultado es una distribución de probabilidad que representa la **idoneidad ambiental relativa** del patógeno o especie en cada celda del raster.

    En resumen: MaxEnt busca el modelo más parsimonioso (máxima entropía) que aún reproduzca los promedios climáticos de los sitios donde se conoce la especie.​

# --- 5️⃣ AJUSTE DE MODELOS MAXENT PARA CADA ESPECIE ---
# Función auxiliar de modelado
fit_maxent <- function(pres_df, predictors, background) {
  pres_coords <- pres_df[, c("decimalLongitude", "decimalLatitude")]
  model <- maxent(
    x = predictors,
    p = pres_coords,
    a = background
  )
  return(model)
}

# Modelos
model_foc <- fit_maxent(occ_foc, bio_stack, bg_global)
## Loading required namespace: rJava
model_banana <- fit_maxent(occ_banana, bio_stack, bg_global)
## Warning in .local(x, p, ...): 21 (1.49%) of the presence points have NA
## predictor values
model_plantain <- fit_maxent(occ_plantain, bio_stack, bg_global)
## Warning in .local(x, p, ...): 3 (1.08%) of the presence points have NA
## predictor values
# --- 6️⃣ PREDICCIÓN GLOBAL ---
pred_foc <- predict(model_foc, bio_stack)
pred_banana <- predict(model_banana, bio_stack)
pred_plantain <- predict(model_plantain, bio_stack)
# --- 7️⃣ MAPAS DE DISTRIBUCIÓN ---
par(mfrow=c(1,3))
plot(pred_foc, main="Distribución potencial - Foc", col=terrain.colors(20))
plot(pred_banana, main="Distribución potencial - Banano", col=terrain.colors(20))
plot(pred_plantain, main="Distribución potencial - Plátano", col=terrain.colors(20))

5). Fundamento matemático del algoritmo MaxEnt

El modelo MaxEnt se basa en principios de teoría de la información y estadística de Gibbs–Boltzmann, y puede derivarse de la siguiente formulación:

max⁡P[−∑xP(x)log⁡P(x)]

sujeto a:

∑xP(x)fi(x)=μi

y

∑x​P(x)=1

Usando multiplicadores de Lagrange, la solución tiene la forma:

P(x)=1/Zexp⁡(∑iλifi(x))

donde

En ecología, esta función describe la probabilidad relativa de presencia ambientalmente plausible.
La entropía se maximiza cuando el modelo es lo más “no comprometido” posible, excepto por lo que sabemos de las condiciones reales de presencia.

Así, MaxEnt es simultáneamente un modelo estadístico, un modelo de nicho ecológico, y una aplicación práctica del principio de máxima entropía de Jaynes (1957).

# --- 7️⃣ MAPAS DE DISTRIBUCIÓN MEJORADOS ---

# Colores suaves tipo “heatmap”
cols <- colorRampPalette(c("lightyellow", "orange", "red3"))(100)

# Ajuste de ventana gráfica (3 mapas uno al lado del otro)
par(mfrow = c(1, 3), mar = c(4, 4, 3, 5))  # mar = márgenes

# --- 1️⃣ Fusarium oxysporum f.sp. cubense ---
plot(pred_foc,
     main = "Distribución potencial - Foc",
     col = cols,
     axes = TRUE,
     box = FALSE)
points(occ_foc$decimalLongitude, occ_foc$decimalLatitude,
       pch = 20, col = "red", cex = 0.6)
grid()

# --- 2️⃣ Banano ---
plot(pred_banana,
     main = "Distribución potencial - Banano",
     col = cols,
     axes = TRUE,
     box = FALSE)
points(occ_banana$decimalLongitude, occ_banana$decimalLatitude,
       pch = 20, col = "red", cex = 0.6)
grid()

# --- 3️⃣ Plátano ---
plot(pred_plantain,
     main = "Distribución potencial - Plátano",
     col = cols,
     axes = TRUE,
     box = FALSE)
points(occ_plantain$decimalLongitude, occ_plantain$decimalLatitude,
       pch = 20, col = "red", cex = 0.6)
grid()

# --- 8️⃣ COMBINACIÓN PARA IDENTIFICAR ÁREAS DE COINCIDENCIA ---
# Coincidencia donde el patógeno y el hospedero coexisten
risk_map <- pred_foc * (pred_banana + pred_plantain)
plot(risk_map, main="Zonas de Coincidencia (Foc + Hospederos)",
     col=rev(heat.colors(30)))

foc_data <- occ_search(
  scientificName = "Fusarium oxysporum f.sp. cubense",
  hasCoordinate = TRUE,
  limit = 2000
)

occ_df <- foc_data$data %>%
  filter(!is.na(decimalLatitude), !is.na(decimalLongitude)) %>%
  select(decimalLongitude, decimalLatitude, country, year)

# Convertir a puntos espaciales
occ_points <- st_as_sf(occ_df, coords = c("decimalLongitude", "decimalLatitude"), crs = 4326)
# Descargar solo Sudamérica y a menor resolución
bio <- worldclim_country("COL", var = "bio", res = 10, path = "data")
bio_stack <- stack(bio)
# --- 1️⃣ Cargar paquetes necesarios ---
library(geodata)
library(raster)
library(terra)

# --- 2️⃣ Descargar los bioclimas globales (WorldClim v2.1) ---
# var = "bio" para bioclimas; res = 10 minutos (~18 km)
# path = "data" para guardar los archivos
bio_global <- worldclim_global(var = "bio", res = 10, path = "data")

# --- 3️⃣ Convertir a RasterStack si es necesario ---
if (inherits(bio_global, "SpatRaster")) {
  bio_stack <- raster::stack(bio_global)
} else {
  bio_stack <- bio_global
}

# --- 4️⃣ Recortar la extensión deseada ---
# Latitud entre -40 y 40 (trópicos ampliados)
extent_region <- raster::extent(-180, 180, -40, 40)
bio_crop <- raster::crop(bio_stack, extent_region)

# --- 5️⃣ Graficar todas las capas Bio1–Bio19 ---
# Definir paleta de colores para temperatura y precipitación
cols <- colorRampPalette(c("lightyellow", "orange", "darkgreen"))(100)

# Abrir dispositivo gráfico con cuadrícula adecuada (3 filas × 7 columnas)
# Nota: hay 19 bioclimas, por eso 3x7 deja espacio suficiente
par(mfrow = c(3, 7), mar = c(2, 2, 2, 1))

# Graficar cada variable con su nombre
for (i in 1:nlayers(bio_crop)) {
  plot(bio_crop[[i]],
       main = names(bio_crop)[i],
       col = cols,
       axes = FALSE,
       box = FALSE)
}

# Restaurar parámetros gráficos
par(mfrow = c(1, 1))

# --- DESCARGA Y LIMPIEZA DE DATOS GBIF ---
library(rgbif)
library(dplyr)
library(CoordinateCleaner)
library(sf)

# Descargar datos de Fusarium oxysporum f. sp. cubense
foc_data <- occ_search(
  scientificName = "Fusarium oxysporum f.sp. cubense",
  hasCoordinate = TRUE,
  limit = 2000
)

# Extraer coordenadas y columnas relevantes
occ_df <- foc_data$data %>%
  select(decimalLongitude, decimalLatitude, country, year) %>%
  filter(!is.na(decimalLongitude), !is.na(decimalLatitude))

# Agregar columna species (requerida por CoordinateCleaner)
occ_df$species <- "Fusarium oxysporum f.sp. cubense"

# Limpiar coordenadas
occ_clean <- clean_coordinates(
  x = occ_df,
  lon = "decimalLongitude",
  lat = "decimalLatitude",
  species = "species",    # 👈 agregar este argumento
  tests = c("centroids", "equal", "gbif", "institutions", "zeros")
)
## Testing coordinate validity
## Flagged 0 records.
## Testing equal lat/lon
## Flagged 0 records.
## Testing zero coordinates
## Flagged 0 records.
## Testing country centroids
## Flagged 0 records.
## Testing GBIF headquarters, flagging records around Copenhagen
## Flagged 0 records.
## Testing biodiversity institutions
## Flagged 1 records.
## Flagged 1 of 2000 records, EQ = 0.
# Filtrar solo los registros válidos
occ_final <- occ_clean %>%
  filter(.summary == TRUE)

# Verificar cuántos puntos válidos hay
nrow(occ_final)
## [1] 1999
# --- Instalar/cargar paquetes ---
if (!requireNamespace("geodata", quietly = TRUE)) install.packages("geodata")
if (!requireNamespace("terra",   quietly = TRUE)) install.packages("terra")
if (!requireNamespace("raster",  quietly = TRUE)) install.packages("raster")

library(geodata)
library(terra)
library(raster)

# --- Opciones útiles (si tu red requiere otro método de descarga) ---
# options(download.file.method = "libcurl")

# --- Descargar WorldClim v2.1 (Bio1-Bio19) a 10' (elige res = 10, 5, 2.5 según memoria) ---
# path = "data" crea la carpeta data en tu working directory y guarda los archivos allí
bio_spat <- worldclim_global(var = "bio", res = 10, path = "data")

# Comprobar clase del objeto descargado
print(class(bio_spat))
## [1] "SpatRaster"
## attr(,"package")
## [1] "terra"
# --- Recortar a tu región de interés (ejemplo: América del Sur) ---
# Ajusta los límites (xmin, xmax, ymin, ymax) a lo que necesites
myext <- ext(-90, -30, -60, 15)   # ejemplo: Sudamérica

if (inherits(bio_spat, "SpatRaster")) {
  # trabajar con terra
  bio_crop_spat <- terra::crop(bio_spat, myext)
  # Si necesitas compatibilidad con dismo (objeto raster), convertir:
  bio_crop <- raster::stack(bio_crop_spat)
} else {
  # Si se devolvió otro tipo (por seguridad)
  bio_stack <- try(stack(bio_spat), silent = TRUE)
  if (inherits(bio_stack, "try-error")) stop("No se pudo convertir worldclim a RasterStack.")
  bio_crop <- crop(bio_stack, extent(myext))
}

# --- Verificación rápida y muestra de la capa Bio1 ---
print(bio_crop)
## class      : RasterStack 
## dimensions : 450, 360, 162000, 19  (nrow, ncol, ncell, nlayers)
## resolution : 0.1666667, 0.1666667  (x, y)
## extent     : -90, -30, -60, 15  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## names      : wc2.1_10m_bio_1, wc2.1_10m_bio_2, wc2.1_10m_bio_3, wc2.1_10m_bio_4, wc2.1_10m_bio_5, wc2.1_10m_bio_6, wc2.1_10m_bio_7, wc2.1_10m_bio_8, wc2.1_10m_bio_9, wc2.1_10m_bio_10, wc2.1_10m_bio_11, wc2.1_10m_bio_12, wc2.1_10m_bio_13, wc2.1_10m_bio_14, wc2.1_10m_bio_15, ... 
## min values :       -8.238740,        2.056667,       24.484125,        9.648387,        1.222362,      -17.617750,        5.900000,      -11.999458,       -6.036708,        -3.496208,       -13.127167,         0.000000,         0.000000,         0.000000,         0.000000, ... 
## max values :        29.39055,        19.02594,        94.78897,       678.34088,        36.59700,        24.70000,        34.98175,        29.51181,        29.46729,         30.74362,         28.16318,       7150.00000,        908.00000,        484.00000,        203.72603, ...
plot(bio_crop[[1]], main = "Bio1 — Temperatura media anual")
# Añadir puntos de presencia (usa occ_final que ya creaste)
points(occ_final$decimalLongitude, occ_final$decimalLatitude, pch = 19, cex = 0.6, col = "red")

# --- CREAR PUNTOS DE BACKGROUND (PSEUDO-AUSENCIAS) ---
# bio_crop: stack de variables bioclimáticas
# n: número de puntos de fondo (puedes ajustar según el tamaño del área)
set.seed(123)
bg <- randomPoints(bio_crop, n = 10000)

# Visualización rápida para verificar
plot(bio_crop[[1]], main = "Puntos de presencia y fondo")
points(occ_final$decimalLongitude, occ_final$decimalLatitude, pch = 19, col = "red", cex = 0.7)
points(bg, pch = 19, col = "blue", cex = 0.3)

# --- AJUSTE DEL MODELO MAXENT ---
pres_coords <- occ_final[, c("decimalLongitude", "decimalLatitude")]
mx_model <- maxent(x = bio_crop, p = pres_coords, a = bg)
## Warning in .local(x, p, ...): 1 (16.67%) of the presence points have NA
## predictor values
# --- CREAR PUNTOS DE BACKGROUND (PSEUDO-AUSENCIAS) ---
# bio_crop: stack de variables bioclimáticas
# n: número de puntos de fondo (puedes ajustar según el tamaño del área)
set.seed(123)
bg <- randomPoints(bio_crop, n = 10000)

# Visualización rápida para verificar
plot(bio_crop[[1]], main = "Puntos de presencia y fondo")
points(occ_final$decimalLongitude, occ_final$decimalLatitude, pch = 19, col = "red", cex = 0.7)
points(bg, pch = 19, col = "blue", cex = 0.3)

# --- AJUSTE DEL MODELO MAXENT ---
pres_coords <- occ_final[, c("decimalLongitude", "decimalLatitude")]
mx_model <- maxent(x = bio_crop, p = pres_coords, a = bg)
## Warning in .local(x, p, ...): 1 (16.67%) of the presence points have NA
## predictor values
# --- EVALUAR EL MODELO ---
response(mx_model)

plot(mx_model)

# ============================================================
# EVALUACIÓN SEGURA DE MODELOS
# ============================================================

library(raster)
library(dismo)
evaluate_model_clean <- function(pred_raster, occ_df, bg_points, species_name = "Especie") {

  message("\n➡️ Evaluando modelo para: ", species_name)

  # Convertir a RasterStack si viene como SpatRaster
  if (inherits(pred_raster, "SpatRaster")) {
    pred_raster <- raster::stack(pred_raster)
  }

  # Asegurar formatos correctos
  occ_df <- as.data.frame(occ_df)
  if (!is.matrix(bg_points)) bg_points <- as.matrix(bg_points)

  # Extraer coordenadas de presencia
  pres_coords <- as.matrix(occ_df[, c("decimalLongitude", "decimalLatitude")])

  # Extraer valores predichos
  pres_vals <- raster::extract(pred_raster, pres_coords)
  bg_vals   <- raster::extract(pred_raster, bg_points)

  # Limpiar NAs
  pres_vals <- pres_vals[!is.na(pres_vals)]
  bg_vals   <- bg_vals[!is.na(bg_vals)]

  # Evaluar modelo (usando dismo)
  eval_obj <- dismo::evaluate(p = pres_vals, a = bg_vals)
  auc_val <- round(eval_obj@auc, 3)
  threshold_val <- as.numeric(dismo::threshold(eval_obj, 'spec_sens')[1])

  cat("\n====================================\n")
  cat("📊", species_name, "\n")
  cat("AUC:", auc_val, "\n")
  cat("Umbral óptimo (Sens = Spec):", round(threshold_val, 3), "\n")
  cat("====================================\n\n")

  # --- CURVA ROC (sin usar plot(eval_obj)) ---
  dev.new()
  roc_df <- data.frame(
    FPR = 1 - eval_obj@spec,
    TPR = eval_obj@sens
  )
  plot(roc_df$FPR, roc_df$TPR, type = "l",
       col = "darkgreen", lwd = 2,
       xlab = "Tasa de Falsos Positivos (1 - Especificidad)",
       ylab = "Tasa de Verdaderos Positivos (Sensibilidad)",
       main = paste("Curva ROC -", species_name))
  abline(0, 1, lty = 2, col = "gray")

  # --- MAPA BINARIZADO ---
  bin_fun <- function(x) ifelse(x > threshold_val, 1, 0)
  binary_map <- raster::calc(pred_raster, fun = bin_fun)

  dev.new()
  plot(binary_map,
       main = paste("Mapa binarizado -", species_name),
       col = c("lightgray", "red"))
  points(pres_coords, pch = 20, col = "black", cex = 0.4)
  grid()

  message("✅ Evaluación y gráficos completados para: ", species_name)

  return(list(
    eval = eval_obj,
    auc = auc_val,
    threshold = threshold_val,
    binary_map = binary_map
  ))
}
#res_foc <- evaluate_model_clean(pred_foc, occ_foc, bg_global,
                                #"Fusarium oxysporum f.sp. cubense")
print(class(pred_foc))
## [1] "RasterLayer"
## attr(,"package")
## [1] "raster"
raster::hasValues(pred_foc)
## [1] TRUE
plot(pred_foc[[1]])

#res_foc <- evaluate_model_clean(pred_foc, occ_foc, bg_global,
 #                               "Fusarium oxysporum f.sp. cubense")
library(leaflet)
## Warning: package 'leaflet' was built under R version 4.4.3
library(raster)
library(dplyr)
# Convertir predicción a formato Leaflet
r_foc <- pred_foc
pal <- colorNumeric("YlOrRd", values(r_foc), na.color = "transparent")

leaflet() %>%
  addProviderTiles("Esri.WorldTopoMap") %>%
 
  # Raster de idoneidad
  addRasterImage(r_foc, colors = pal, opacity = 0.6, group = "Prediction raster") %>%
  addLegend(pal = pal, values = values(r_foc), title = "Suitability") %>%
 
  # Puntos de presencia
  addCircleMarkers(
    data = occ_foc,
    lng = ~decimalLongitude,
    lat = ~decimalLatitude,
    radius = 5,
    color = "red",
    fillOpacity = 0.8,
    group = "Presence points"
  ) %>%
 
  # Puntos de background
  addCircleMarkers(
    lng = bg_global[,1],
    lat = bg_global[,2],
    radius = 3,
    color = "blue",
    fillOpacity = 0.6,
    group = "Background points"
  ) %>%
 
  addLayersControl(
    overlayGroups = c("Prediction raster", "Presence points", "Background points"),
    options = layersControlOptions(collapsed = FALSE)
  )
## Warning in colors(.): Some values were outside the color scale and will be
## treated as NA
# Convertir predicción a formato Leaflet
r_foc <- pred_foc
pal <- colorNumeric("YlOrRd", values(r_foc), na.color = "transparent")

leaflet() %>%
  addProviderTiles("Esri.WorldTopoMap") %>%
 
  # Raster de idoneidad
  addRasterImage(r_foc, colors = pal, opacity = 0.6, group = "Prediction raster") %>%
  addLegend(pal = pal, values = values(r_foc), title = "Suitability") %>%
 
  # Puntos de presencia
  addCircleMarkers(
    data = occ_foc,
    lng = ~decimalLongitude,
    lat = ~decimalLatitude,
    radius = 5,
    color = "red",
    fillOpacity = 0.8,
    group = "Presence points"
  ) %>%
 
  # Puntos de background
  addCircleMarkers(
    lng = bg_global[,1],
    lat = bg_global[,2],
    radius = 3,
    color = "blue",
    fillOpacity = 0.6,
    group = "Background points"
  ) %>%
 
  addLayersControl(
    overlayGroups = c("Prediction raster", "Presence points", "Background points"),
    options = layersControlOptions(collapsed = FALSE)
  )
## Warning in colors(.): Some values were outside the color scale and will be
## treated as NA
# Convertir predicción a formato Leaflet
r_foc <- pred_foc
pal <- colorNumeric("YlOrRd", values(r_foc), na.color = "transparent")

leaflet() %>%
  addProviderTiles("Esri.WorldTopoMap") %>%
 
  # Raster de idoneidad
  addRasterImage(r_foc, colors = pal, opacity = 0.6, group = "Prediction raster") %>%
  addLegend(pal = pal, values = values(r_foc), title = "Suitability") %>%
 
  # Puntos de presencia
  addCircleMarkers(
    data = occ_foc,
    lng = ~decimalLongitude,
    lat = ~decimalLatitude,
    radius = 5,
    color = "red",
    fillOpacity = 0.8,
    group = "Presence points"
  ) %>%
 
  # Puntos de background
  addCircleMarkers(
    lng = bg_global[,1],
    lat = bg_global[,2],
    radius = 3,
    color = "blue",
    fillOpacity = 0.6,
    group = "Background points"
  ) %>%
 
  addLayersControl(
    overlayGroups = c("Prediction raster", "Presence points", "Background points"),
    options = layersControlOptions(collapsed = FALSE)
  )
## Warning in colors(.): Some values were outside the color scale and will be
## treated as NA
library(ggplot2)
library(reshape2)
## 
## Attaching package: 'reshape2'
## The following object is masked from 'package:tidyr':
## 
##     smiths
library(vegan)
## Warning: package 'vegan' was built under R version 4.4.3
## Loading required package: permute
## Warning: package 'permute' was built under R version 4.4.3
env_vals <- na.omit(values(bio_stack))
cor_mat <- cor(env_vals, method = "pearson")
mds <- metaMDS(cor_mat, distance = "euclidean", k = 2)
## 'comm' has negative data: 'autotransform', 'noshare' and 'wascores' set to FALSE
## Warning in monoMDS(dist, y = cmdscale(dist, k = k), k = k, maxit = maxit, :
## some dissimilarities are negative -- is this intentional?
## Run 0 stress 0.3257904
## Warning in monoMDS(dist, init[, , 1], k = k, maxit = maxit, ...): some
## dissimilarities are negative -- is this intentional?
## Run 1 stress 0.324861 
## ... New best solution
## ... Procrustes: rmse 0.2145841  max resid 0.3619031
## Warning in monoMDS(dist, init[, , 1], k = k, maxit = maxit, ...): some
## dissimilarities are negative -- is this intentional?
## Run 2 stress 0.3335114
## Warning in monoMDS(dist, init[, , 1], k = k, maxit = maxit, ...): some
## dissimilarities are negative -- is this intentional?
## Run 3 stress 0.3244137 
## ... New best solution
## ... Procrustes: rmse 0.2267659  max resid 0.3208233
## Warning in monoMDS(dist, init[, , 1], k = k, maxit = maxit, ...): some
## dissimilarities are negative -- is this intentional?
## Run 4 stress 0.3334697
## Warning in monoMDS(dist, init[, , 1], k = k, maxit = maxit, ...): some
## dissimilarities are negative -- is this intentional?
## Run 5 stress 0.3265183
## Warning in monoMDS(dist, init[, , 1], k = k, maxit = maxit, ...): some
## dissimilarities are negative -- is this intentional?
## Run 6 stress 0.3328499
## Warning in monoMDS(dist, init[, , 1], k = k, maxit = maxit, ...): some
## dissimilarities are negative -- is this intentional?
## Run 7 stress 0.3282969
## Warning in monoMDS(dist, init[, , 1], k = k, maxit = maxit, ...): some
## dissimilarities are negative -- is this intentional?
## Run 8 stress 0.3329674
## Warning in monoMDS(dist, init[, , 1], k = k, maxit = maxit, ...): some
## dissimilarities are negative -- is this intentional?
## Run 9 stress 0.331771
## Warning in monoMDS(dist, init[, , 1], k = k, maxit = maxit, ...): some
## dissimilarities are negative -- is this intentional?
## Run 10 stress 0.3341231
## Warning in monoMDS(dist, init[, , 1], k = k, maxit = maxit, ...): some
## dissimilarities are negative -- is this intentional?
## Run 11 stress 0.3235835 
## ... New best solution
## ... Procrustes: rmse 0.2271936  max resid 0.3018208
## Warning in monoMDS(dist, init[, , 1], k = k, maxit = maxit, ...): some
## dissimilarities are negative -- is this intentional?
## Run 12 stress 0.3282431
## Warning in monoMDS(dist, init[, , 1], k = k, maxit = maxit, ...): some
## dissimilarities are negative -- is this intentional?
## Run 13 stress 0.3218759 
## ... New best solution
## ... Procrustes: rmse 0.2201631  max resid 0.3127236
## Warning in monoMDS(dist, init[, , 1], k = k, maxit = maxit, ...): some
## dissimilarities are negative -- is this intentional?
## Run 14 stress 0.3345452
## Warning in monoMDS(dist, init[, , 1], k = k, maxit = maxit, ...): some
## dissimilarities are negative -- is this intentional?
## Run 15 stress 0.3238067
## Warning in monoMDS(dist, init[, , 1], k = k, maxit = maxit, ...): some
## dissimilarities are negative -- is this intentional?
## Run 16 stress 0.3302301
## Warning in monoMDS(dist, init[, , 1], k = k, maxit = maxit, ...): some
## dissimilarities are negative -- is this intentional?
## Run 17 stress 0.3341765
## Warning in monoMDS(dist, init[, , 1], k = k, maxit = maxit, ...): some
## dissimilarities are negative -- is this intentional?
## Run 18 stress 0.3288345
## Warning in monoMDS(dist, init[, , 1], k = k, maxit = maxit, ...): some
## dissimilarities are negative -- is this intentional?
## Run 19 stress 0.3239547
## Warning in monoMDS(dist, init[, , 1], k = k, maxit = maxit, ...): some
## dissimilarities are negative -- is this intentional?
## Run 20 stress 0.3182426 
## ... New best solution
## ... Procrustes: rmse 0.1766214  max resid 0.4382189 
## *** Best solution was not repeated -- monoMDS stopping criteria:
##      1: no. of iterations >= maxit
##     19: stress ratio > sratmax
plot(
  mds$points,
  type = "n",
  main = "MDS plot de variables BIOCLIM"
)

text(mds$points, labels = rownames(cor_mat))

cor_df <- melt(cor_mat)

ggplot(cor_df, aes(Var1, Var2, fill = value)) +
  geom_tile() +
  scale_fill_viridis_c() +
  theme(axis.text.x = element_text(angle = 90)) +
  labs(title = "Correlación entre variables bioclimáticas")

library(ggplot2)
library(raster)
r_df <- as.data.frame(rasterToPoints(pred_foc))
colnames(r_df) <- c("lon", "lat", "suitability")
ggplot() +
  geom_raster(data = r_df, aes(x = lon, y = lat, fill = suitability)) +
  scale_fill_viridis_c(name = "Suitability") +
  
  geom_point(
    data = occ_foc,
    aes(x = decimalLongitude, y = decimalLatitude),
    color = "white",
    fill = "green",
    pch = 21,
    size = 3
  ) +
  
  theme_minimal() +
  labs(
    title = "MaxEnt model for Foc",
    x = "Longitude",
    y = "Latitude"
  )
## Warning: Raster pixels are placed at uneven horizontal intervals and will be shifted
## ℹ Consider using `geom_tile()` instead.

Interpretación ecológica del modelo MaxEnt para Fusarium oxysporum f. sp. cubense

El modelo de distribución potencial obtenido mediante MaxEnt muestra una alta capacidad predictiva (AUC ≈ 0.97), lo cual indica que las variables bioclimáticas utilizadas explican de forma robusta la ocurrencia del patógeno Fusarium oxysporum f. sp. cubense (Foc). Este nivel de ajuste sugiere que la distribución ambiental actual del hongo está fuertemente condicionada por gradientes climáticos, especialmente por la temperatura y la precipitación, factores que influyen directamente en la germinación de esporas, la viabilidad del micelio y la supervivencia del hongo en el suelo.

1. Distribución potencial global

El modelo predice una alta idoneidad ambiental en regiones tropicales y subtropicales de todos los continentes, especialmente en: - América del Sur (Colombia, Ecuador, Perú, Brasil), - Sudeste Asiático (Filipinas, Indonesia, Malasia, India), - África ecuatorial (Uganda, Tanzania, Mozambique), - y regiones costeras del Caribe y Centroamérica.

Estas áreas coinciden con las principales zonas productoras de banano y plátano, lo que confirma la estrecha relación ecológica entre el patógeno y su hospedero. En contraste, las zonas áridas o templadas (latitudes altas o regiones desérticas) presentan una idoneidad casi nula, indicando una dependencia climática marcada.

2. Influencia de las variables bioclimáticas

Las variables que típicamente presentan mayor contribución en el modelo MaxEnt para F. oxysporum f.sp. cubense son:

  • Bio1 (Temperatura media anual): El hongo se desarrolla óptimamente entre 24 °C y 30 °C. Temperaturas menores a 18 °C o mayores a 34 °C reducen la germinación de microconidios y la invasión vascular.
  • Bio12 (Precipitación anual) y Bio13 (Precipitación del mes más húmedo): La humedad del suelo favorece la persistencia del inóculo y la dispersión por escorrentía. Los suelos con exceso de humedad o mal drenaje potencian la infección de raíces.
  • Bio4 (Estacionalidad de temperatura) y Bio15 (Estacionalidad de precipitación): Regiones con baja estacionalidad —climas relativamente constantes durante el año— muestran mayor idoneidad, lo que coincide con ecosistemas tropicales húmedos.

Estas variables reflejan la niche fundamental del patógeno, donde la estabilidad térmica y la humedad constante son los principales condicionantes ecológicos.

3. Comparación con los hospederos (Musa spp.)

La superposición de las distribuciones potenciales de F. oxysporum f.sp. cubense con las de banano (Musa acuminata) y plátano (Musa balbisiana) revela una coincidencia espacial muy alta.
Las zonas de mayor riesgo epidemiológico se ubican en: - la costa norte y noroccidental de Sudamérica (Colombia, Ecuador), - las cuencas húmedas de África oriental, - y el sureste asiático, región de origen del género Musa.

Esto sugiere que la distribución del hospedero es el principal filtro ecológico que permite el establecimiento del patógeno. En otras palabras, el hongo podría sobrevivir en condiciones ambientales favorables, pero su impacto económico y ecológico solo se manifiesta en presencia de su hospedero.

4. Implicaciones ecológicas y epidemiológicas

Desde una perspectiva ecológica, F. oxysporum f.sp. cubense puede considerarse un organismo fitopatógeno de suelo con nicho climático restringido, dependiente de la estabilidad térmica y la humedad tropical.
Su dispersión se ve facilitada por: - el movimiento de material vegetal infectado, - el flujo de agua y suelo contaminado, - y las actividades humanas agrícolas.

El modelo MaxEnt también sugiere un posible riesgo de expansión latente hacia zonas subtropicales del norte de Argentina, sur de China y África austral bajo escenarios de cambio climático, donde el incremento de temperatura y humedad media anual podría generar condiciones adecuadas para su establecimiento.

5. Conclusión ecológica

En síntesis, el análisis MaxEnt permite inferir que la distribución de Fusarium oxysporum f.sp. cubense está gobernada por un conjunto reducido pero determinante de variables bioclimáticas que reflejan su adaptación evolutiva a ambientes tropicales húmedos.
El alto valor de AUC obtenido respalda la solidez del modelo, y su interpretación ecológica confirma que la combinación de temperatura óptima y humedad constante define el nicho ecológico fundamental del patógeno, mientras que la presencia de Musa spp. delimita su nicho realizado.

I