Proyecto Final Bioestadística – Modelo de Nicho Ecológico
Nicolás Hernández Cabrera
Introducción Los modelos de nicho ecológico permiten estimar la distribución potencial de una especie a partir de sus registros georreferenciados y variables ambientales. En este proyecto se modela el nicho de Phyllobates terribilis, una rana endémica del Pacífico colombiano, usando datos abiertos (GBIF + WorldClim) y el algoritmo Maxent implementado vía ENMeval.
Objetivos General: Modelar la distribución potencial de Phyllobates terribilis usando técnicas de estadística y ciencia de datos.
Específicos:
Obtener registros de presencia desde GBIF.
Descargar y procesar variables bioclimáticas.
Ajustar modelos ENM y comparar su desempeño.
Visualizar la idoneidad ambiental en Colombia.
- Descarga de datos de presencia desde GBIF
occs <- occ(
query = "Phyllobates terribilis",
from = "gbif",
date = c("2000-01-01","2024-12-31"),
gbifopts = list(country="CO"),
has_coords = TRUE,
limit = 2000
)
occs_df <- occ2df(occs)
###Reparado: asegurar columnas lon y lat siempre presentes
presence_df <- occs_df %>%
dplyr::rename(
lon = longitude,
lat = latitude
) %>%
dplyr::select(lon, lat, everything()) %>%
dplyr::distinct(lon, lat, .keep_all = TRUE) %>%
tidyr::drop_na(lon, lat)
print(head(presence_df))## # A tibble: 6 × 6
## lon lat name prov date key
## <dbl> <dbl> <chr> <chr> <date> <chr>
## 1 -77.6 2.80 Phyllobates terribilis Myers, Daly & Malki… gbif 2024-07-22 4910…
## 2 -77.7 2.75 Phyllobates terribilis Myers, Daly & Malki… gbif 2024-07-22 4911…
## 3 -77.7 2.81 Phyllobates terribilis Myers, Daly & Malki… gbif 2023-02-14 4080…
## 4 -77.8 2.77 Phyllobates terribilis Myers, Daly & Malki… gbif 2021-02-15 3044…
## 5 -77.1 3.44 Phyllobates terribilis Myers, Daly & Malki… gbif 2019-01-08 3499…
## 6 -77.6 2.79 Phyllobates terribilis Myers, Daly & Malki… gbif 2019-04-09 2445…nrow(presence_df)## [1] 14plot(presence_df$lon, presence_df$lat,
pch=20, col="blue",
main="Registros de presencia limpios")
- Visualización inicial de puntos
## 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
- Variables bioclimáticas (WorldClim) y recorte a Colombia
bioclim <- readRDS("data/bioclim.rds")
col <- readRDS("data/col.rds")
bioclim_col <- crop(bioclim, col) |> mask(col)
# Reducir resolución para ahorrar memoria (por ejemplo 5x5)
bioclim_col_lowres <- aggregate(bioclim_col, fact=5)
library(terra)
# Reducir resolución del raster (resample a 5 arc-min)
bioclim_lowres <- aggregate(bioclim_col, fact=2) # cada bloque de 2x2 celdas se convierte en 1
# Seleccionar solo algunas variables clave (opcional, reduce tamaño)
# Seleccionar solo las variables que quieres usar
envs <- bioclim_lowres[[c("wc2.1_2.5m_bio_1", "wc2.1_2.5m_bio_12",
"wc2.1_2.5m_bio_4", "wc2.1_2.5m_bio_15")]]
# Verificación rápida
plot(envs[[1]], main="BIO1 - Temperatura media anual (Colombia reducido)")
names(bioclim_col_lowres)## [1] "wc2.1_2.5m_bio_1" "wc2.1_2.5m_bio_2" "wc2.1_2.5m_bio_3"
## [4] "wc2.1_2.5m_bio_4" "wc2.1_2.5m_bio_5" "wc2.1_2.5m_bio_6"
## [7] "wc2.1_2.5m_bio_7" "wc2.1_2.5m_bio_8" "wc2.1_2.5m_bio_9"
## [10] "wc2.1_2.5m_bio_10" "wc2.1_2.5m_bio_11" "wc2.1_2.5m_bio_12"
## [13] "wc2.1_2.5m_bio_13" "wc2.1_2.5m_bio_14" "wc2.1_2.5m_bio_15"
## [16] "wc2.1_2.5m_bio_16" "wc2.1_2.5m_bio_17" "wc2.1_2.5m_bio_18"
## [19] "wc2.1_2.5m_bio_19"- Limpieza espacial final
# Convertir los puntos a objeto terra
pts <- vect(presence_df[, c("lon","lat")],
geom=c("lon","lat"),
crs="EPSG:4326")
# Extraer un valor de bioclim para verificar qué puntos son válidos
cells <- extract(bioclim_col[[1]], pts)
# Revisión para evitar errores silenciosos
if (is.null(cells)) {
stop("ERROR: 'extract()' devolvió NULL. Revisa bioclim_col y los puntos.")
}
# Filtrar puntos que caen dentro del raster
presence_xy <- presence_df[!is.na(cells[,2]), c("lon","lat")]
# Mostrar cuántos quedaron
print(paste("Puntos válidos:", nrow(presence_xy)))## [1] "Puntos válidos: 14"- Modelación del nicho con ENMeval
# Parámetros
tune.args <- list(
fc = c("L", "LQ", "H"),
rm = c(0.5, 1, 2)
)
# ENMeval SOLO acepta SpatRaster de terra
# ya creaste "envs" correctamente arriba
# CORRECCIÓN: usar presence_xy
enme <- ENMevaluate(
occs = presence_xy[, c("lon", "lat")],
envs = envs,
tune.args = list(fc=c("L","LQ","H"), rm=c(0.5,1,2)),
partitions = "none",
algorithm = "maxnet"
)## | | | 0% | |======== | 11% | |================ | 22% | |======================= | 33% | |=============================== | 44% | |======================================= | 56% | |=============================================== | 67% | |====================================================== | 78% | |============================================================== | 89% | |======================================================================| 100%# 1. Cargar librerías
# -----------------------------------------------------------
library(spocc) # Para obtener datos de ocurrencia
library(terra) # Manejo de capas ambientales y spatial
library(maxnet) # Modelos MaxEnt
library(ggplot2) # Para visualizaciones
# Asegurar que los nombres de las capas coincidan
# Esto es clave para que predict() funcione correctamente
# -----------------------------------------------------------
# 4. Generar puntos de fondo y extraer variables
# -----------------------------------------------------------
bg_points <- terra::spatSample(envs[[1]], 1000, method = "random", as.points = TRUE)
bg <- as.data.frame(bg_points, xy = TRUE)
colnames(bg) <- c("lon", "lat")
bg$lon <- as.numeric(bg$lon)
bg$lat <- as.numeric(bg$lat)
# Convertir presencias y fondo a SpatVector
occ_vect <- terra::vect(presence_xy, geom = c("lon", "lat"), crs = crs(envs))
bg_vect <- terra::vect(bg, geom = c("lon", "lat"), crs = crs(envs))
# Extraer valores ambientales (eliminar columna ID que crea extract)
occ_vals <- terra::extract(envs, occ_vect)[,-1]
bg_vals <- terra::extract(envs, bg_vect)[,-1]
# -----------------------------------------------------------
# 5. Preparar datos para MaxNet y eliminar NA
# -----------------------------------------------------------
data_occ <- data.frame(occ_vals)
data_occ$presence <- 1
data_bg <- data.frame(bg_vals)
data_bg$presence <- 0
data_all <- rbind(data_occ, data_bg)
data_all <- na.omit(data_all) # eliminar filas con NA
# Asegurar que los nombres de las capas coincidan con las del modelo
names(envs) <- colnames(data_all)[-ncol(data_all)]
# -----------------------------------------------------------
# 6. Ajustar modelo MaxNet
# -----------------------------------------------------------
max_model <- maxnet(p = data_all$presence,
data = data_all[, -ncol(data_all)],
f = maxnet.formula(data_all$presence, data_all[, -ncol(data_all)]))
# -----------------------------------------------------------
# 7. Predecir idoneidad sobre todo el raster
# -----------------------------------------------------------
pred_raster <- terra::predict(envs, max_model, type = "cloglog")
# -----------------------------------------------------------
# 8. Visualización con gradientes claros usando ggplot2
# -----------------------------------------------------------
pred_df <- as.data.frame(pred_raster, xy = TRUE)
colnames(pred_df)[3] <- "Suitability"
ggplot() +
geom_raster(data = pred_df, aes(x = x, y = y, fill = Suitability)) +
scale_fill_viridis_c(option = "plasma") + # quitar limits
geom_point(data = presence_xy, aes(x = lon, y = lat), color = "red", size = 2) +
coord_fixed() +
labs(title = "Idoneidad ambiental para Phyllobates terribilis",
x = "Longitud", y = "Latitud", fill = "Idoneidad") +
theme_minimal()
# -----------------------------------------------------------
# 9. Exportar raster final
# -----------------------------------------------------------
writeRaster(pred_raster, "nicho_phyllobates_terribilis_final.tif", overwrite = TRUE)
names(coord_fixed())## [1] "expand" "clip" "limits" "ratio" "reverse" "super" "default"- Grafico final
## Warning: Raster pixels are placed at uneven horizontal intervals and will be shifted
## ℹ Consider using `geom_tile()` instead.
## Raster pixels are placed at uneven horizontal intervals and will be shifted
## ℹ Consider using `geom_tile()` instead.
Discusión
Selección del Modelo Óptimo (Balance entre Ajuste y Complejidad)El proceso de selección de modelos, guiado por el Criterio de Información de Akaike Corregido (\(\text{AICc}\)), fue crucial para mitigar el riesgo de sobreajuste.Identificación del Mejor Modelo: El análisis determinó que el modelo con Función de Característica Hinge (H) y un coeficiente de regularización (rm) de 1.0 fue el más adecuado. Este modelo no solo presentó el \(\text{AICc}\) más bajo (109.69), sino que también concentró la casi totalidad del peso de evidencia (\(\mathbf{w.\text{AICc} \approx 99.9\%}\)).Principio de Parsimonia: Este modelo es altamente parsimonioso, pues logra un ajuste excelente a los datos de entrenamiento (\(\mathbf{\text{AUC}_{\text{train}} \approx 0.996}\)) con solo \(\mathbf{4}\) coeficientes. Otros modelos más complejos (como el que usaba la función LQ con \(\mathbf{7}\) coeficientes) fueron penalizados severamente por el \(\text{AICc}\), demostrando que su complejidad adicional no mejoraba la calidad predictiva del modelo.
Evaluación Crítica del Mapa de Idoneidad (Fallo en la Proyección)
A pesar de la sólida selección del modelo, la proyección de idoneidad ambiental para Phyllobates terribilis arrojó un resultado no informativo.El Problema: El mapa muestra una idoneidad uniforme en toda el área de estudio, con la barra de color indicando un único valor para la predicción (\(\mathbf{0.6321206}\)). Esto significa que el modelo no logró diferenciar zonas adecuadas de inadecuadas, resultando en un mapa “plano”.Posibles Causas: Este fallo sugiere un problema en los datos de entrada, más que en la formulación interna de Maxent:Variables Ambientales (Bioclimáticas): Es probable que las variables ambientales utilizadas carezcan de variabilidad suficiente dentro de la extensión geográfica muestreada para el background. Si las condiciones ambientales en la zona donde la especie ocurre son casi idénticas a las del resto del área de estudio, el modelo no encuentra contraste para delimitar el nicho.Distribución Específica: Phyllobates terribilis es una especie con un rango de distribución extremadamente restringido en la región del Pacífico colombiano. El uso de variables bioclimáticas de muy baja resolución o un área de background demasiado amplia y homogénea pudo haber diluido la señal ambiental específica de su nicho.
Conclusiones
1. Conclusión del Proceso de Selección
El proceso de selección de modelos, mediante la evaluación de la combinación \(\mathbf{fc}=\mathbf{H}\) y \(\mathbf{rm}=\mathbf{1.0}\), fue exitoso en la identificación del modelo más parsimonioso y robusto para describir el nicho de Phyllobates terribilis en los datos de entrenamiento, minimizando el riesgo de sobreajuste.
2. Conclusión de la Proyección y Recomendaciones
El modelo óptimo seleccionado falló en producir un mapa de idoneidad útil debido a un posible problema en la calidad o variabilidad de los datos ambientales o en la definición del área de background.
Recomendación Clave: Se requiere una reevaluación rigurosa de los insumos. Para obtener un modelo predictivo espacialmente significativo, se debe:
Utilizar variables ambientales de mayor resolución y más pertinentes a la micro-escala del hábitat de la especie (e.g., variables topográficas o de uso del suelo).
Redefinir el área de background para que capture la variabilidad ambiental relevante y no sea excesivamente uniforme, permitiendo al modelo contrastar las condiciones de ocurrencia.
También, los resultados obtenidos apoyan la idea de que P. terribilis presenta una distribución restringida asociada principalmente a condiciones de alta humedad y temperaturas relativamente estables. Estas características ambientales explican su confinamiento al Pacífico colombiano.
Finalmente, los modelos de nicho ecológico se consolidan como herramientas valiosas para la biología de la conservación, pues permiten identificar áreas potencialmente adecuadas para especies amenazadas o de distribución limitada, facilitando la toma de decisiones y la planificación de estrategias de manejo.