11 Modelos de nichos ecológicos

11.1 Introducción

Los modelos de nichos ecológicos permiten predecir la distribución potencial de especies en función de datos de presencia y variables ambientales (ej. clima, suelo, vegetación, altitud). Puede consultar estas diapositivas sobre el tema.

11.2 Instalación y carga de paquetes

Se utilizan varios paquetes para manejar datos geoespaciales y para ejecutar la modelización.

# Paquete para acceder datos en GBIF
install.packages("rgbif")

# Paquete para acceder datos geoespaciales
install.packages("geodata")

# Paquete para mapas interactivos
install.packages("leaflet")

# Paquete para modelado de distribución de especies
install.packages("dismo")

# Colección de paquetes de Tidyverse
library(tidyverse)

# Estilos para ggplot2
library(ggthemes)

# Paletas de colores de RColorBrewer
library(RColorBrewer)

# Paletas de colores de viridis
library(viridisLite)

# Gráficos interactivos
library(plotly)

# Manejo de datos vectoriales
library(sf)

# Manejo de datos raster
library(terra)

# Manejo de datos raster
library(raster)

# Mapas interactivos
library(leaflet)

# Acceso a datos en GBIF
library(rgbif)

# Datos geoespaciales
library(geodata)

# Modelado de distribución de especies
library(dismo)

11.3 Obtención de datos de presencia

En el siguiente bloque se especifica el nombre de la especie que se va a modelar. Como ejemplo inicial se desarrolla un modelo de distribución de Bradypus variegatus (perezoso de tres dedos).

# Nombre de la especie
especie <- "Bradypus variegatus"

Para obtener los datos de presencia se emplea el paquete rgbif, el cual proporciona acceso a los datos agrupados por el Sistema Mundial de Información en Biodiversidad (GBIF).

# Consulta a GBIF
respuesta <- occ_search(
  scientificName = especie, 
  hasCoordinate = TRUE,
  hasGeospatialIssue = FALSE,
  limit = 10000
)

# Extraer datos de presencia
presencia <- respuesta$data

Opcionalmente, se puede guardar el dataframe con los datos de presencia en un archivo CSV, para no tener que repetir la consulta al API de GBIF.

# Guardar los datos de presencia en un archivo CSV
write_csv(presencia, 'presencia.csv')

El siguiente bloque debe ejecutarse si se desea leer los datos desde el archivo CSV.

# Leer los datos de presencia de un archivo CSV
presencia <- read_csv('presencia.csv')

Para manejarlo como un conjunto de datos geoespacial, el dataframe retornado por occ_search() se convierte a un objeto sf.

presencia <- st_as_sf(
  presencia,
  coords = c("decimalLongitude", "decimalLatitude"),
  remove = FALSE, # conservar las columnas de las coordenadas
  crs = 4326
)

11.3.1 Registros por país

Código

# Gráfico ggplot2
grafico_ggplot2 <-
  presencia |>
  st_drop_geometry() |>
  ggplot(aes(x = fct_infreq(countryCode))) +
  geom_bar(
    aes(
      text = paste0(
        "Cantidad de registros de presencia: ", after_stat(count)
      )
    )    
  ) +
  ggtitle("Cantidad de registros de presencia por país") +
  xlab("País") +
  ylab("Cantidad de registros de presencia") +
  labs(caption = "Fuente: GBIF") +
  theme_economist()

# Gráfico plotly
ggplotly(grafico_ggplot2, tooltip = "text") |> 
  config(locale = 'es')

11.3.2 Registros por año

Código

# Gráfico ggplot2
grafico_ggplot2 <-
  presencia |>
  st_drop_geometry() |>
  group_by(year) |>
  summarize(n = n()) |>
  ggplot(aes(x = year, y = n)) +
  geom_line() +
  geom_point(
    aes(
      text = paste0(
        "Año: ", year, "\n",
        "Cantidad de registros: ", n
      )
    )
  ) +
  ggtitle("Cantidad de registros de presencia por año") +
  xlab("Año") +
  ylab("Cantidad de registros de presencia") +
  labs(caption = "Fuente: GBIF") +
  theme_economist()

# Gráfico plotly
ggplotly(grafico_ggplot2, tooltip = "text") |> 
  config(locale = 'es')

11.3.3 Mapa

Código

# Mapa
leaflet() |>
  addTiles(group = "Mapa general") |>
  addProviderTiles(
    providers$Esri.WorldImagery, 
    group = "Imágenes satelitales"
  ) |>
  addProviderTiles(
    providers$CartoDB.Positron, 
    group = "Mapa blanco"
  ) |>  
  addCircleMarkers(
    # capa de registros de presencia (puntos)
    data = presencia,
    stroke = F,
    radius = 3,
    fillColor = 'red',
    fillOpacity = 1,
    popup = paste(
      paste0("<strong>País: </strong>", presencia$country),
      paste0("<strong>Localidad: </strong>", presencia$locality),
      paste0("<strong>Fecha: </strong>", presencia$eventDate),
      paste0("<strong>Fuente: </strong>", presencia$institutionCode),
      paste0("<a href='", presencia$occurrenceID, "'>Más información</a>"),
      sep = '<br/>'
    ),
    group = "Registros de Bradypus variegatus"
  ) |>
  addLayersControl(
    baseGroups = c("Mapa general", "Imágenes satelitales", "Mapa blanco"),
    overlayGroups = c("Registros de Bradypus variegatus"))

11.4 Obtención de variables ambientales

Frecuentemente se utilizan datos bioclimáticos (ej. temperatura, precipitación) como variables ambientales. Pueden obtenerse mediante el paquete geodata. En este caso, se descargan las variables bioclimáticas de WorldClim.

# Consulta a WorldClim
clima <- worldclim_global(var = 'bio', res = 10, path = tempdir())

# Nombres de las variables climáticas
names(clima)

 [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"

Seguidamente, se “recortan” las capas raster climáticas provenientes de WorldClim para así cubrir solamente el área en la que se encuentra presente la especie.

# Definir la extensión del área de estudio
area_estudio <- ext(
  min(presencia$decimalLongitude) - 5, 
  max(presencia$decimalLongitude) + 5,
  min(presencia$decimalLatitude) - 5, 
  max(presencia$decimalLatitude) + 5
)

# Recortar las variables bioclimáticas al área de estudio
clima <- crop(clima, area_estudio)

11.4.1 Mapa

Código

# Paleta de colores de temperatura
colores_temperatura <- colorNumeric(
  # palette = "inferno",
  # palette = "magma",
  palette = rev(brewer.pal(11, "RdYlBu")),
  values(clima$wc2.1_10m_bio_1),
  na.color = "transparent"
)

# Paleta de colores de precipitación
colores_precipitacion <- colorNumeric(
  # palette = "viridis",
  # palette = "YlGnBu",  
  palette = "Blues",
  values(clima$wc2.1_10m_bio_12),
  na.color = "transparent"
)

# Mapa
leaflet() |>
  addTiles(group = "Mapa general") |>
  addProviderTiles(
    providers$Esri.WorldImagery, 
    group = "Imágenes satelitales"
  ) |>  
  addProviderTiles(
    providers$CartoDB.Positron, 
    group = "Mapa blanco"
  ) |>
  addRasterImage( # capa raster de temperatura
    clima$wc2.1_10m_bio_1,
    colors = colores_temperatura, # paleta de colores
    opacity = 0.6,
    group = "Temperatura",
  ) |>
  addRasterImage( # capa raster de precipitación
    clima$wc2.1_10m_bio_12,
    colors = colores_precipitacion, # paleta de colores
    opacity = 0.6,
    group = "Precipitación",
  ) |>
  addCircleMarkers(
    # capa de registros de presencia (puntos)
    data = presencia,
    stroke = F,
    radius = 3,
    fillColor = 'red',
    fillOpacity = 1,
    popup = paste(
      paste0("<strong>País: </strong>", presencia$country),
      paste0("<strong>Localidad: </strong>", presencia$locality),
      paste0("<strong>Fecha: </strong>", presencia$eventDate),
      paste0("<strong>Fuente: </strong>", presencia$institutionCode),
      paste0("<a href='", presencia$occurrenceID, "'>Más información</a>"),
      sep = '<br/>'
    ),
    group = "Registros de Bradypus variegatus"
  ) |>  
  addLegend(
    title = "Temperatura",
    values = values(clima$wc2.1_10m_bio_1),
    pal = colores_temperatura,
    position = "bottomleft",
    group = "Temperatura"
  ) |>
  addLegend(
    title = "Precipitación",
    values = values(clima$wc2.1_10m_bio_12),
    pal = colores_precipitacion,
    position = "bottomleft",
    group = "Precipitación"
  ) |>  
  addLayersControl(
    # control de capas
    baseGroups = c("Mapa general", "Imágenes satelitales", "Mapa blanco"),
    overlayGroups = c("Temperatura", "Precipitación", "Registros de Bradypus variegatus")
  ) |>
  hideGroup("Precipitación")

11.5 Modelización

11.5.1 Creación de conjuntos de entrenamiento y de evaluación

Primero, se eliminan las coordenadas duplicadas del conjunto de datos de presencia.

# Crear dataframe con columnas de longitud y latitud
coordenadas_presencia <- data.frame(
  decimalLongitude = presencia$decimalLongitude,
  decimalLatitude = presencia$decimalLatitude
)

# Eliminar coordenadas duplicadas
coordenadas_presencia <- unique(coordenadas_presencia)

Seguidamente se dividen los datos de presencia en dos subconjuntos:

Entrenamiento: para desarrollar el modelo.
Evaluación: para evaluar el modelo.

# Establecer una "semilla" para garantizar que la selección aleatoria sea reproducible
set.seed(123)

# Cantidad de registros de presencia
n_presencia <- nrow(coordenadas_presencia)

# Con sample(), se selecciona aleatoriamente una proporción (ej. 0.7) 
# de los índices de los datos de presencia para el conjunto de entrenamiento
indices_entrenamiento <- sample(
  1:n_presencia, 
  size = round(0.7 * n_presencia)
)

# Crear el subconjunto de entrenamiento utilizando los índices seleccionados
entrenamiento <- coordenadas_presencia[indices_entrenamiento, ]

# Crear el subconjunto de evaluación con los datos restantes
evaluacion <- coordenadas_presencia[-indices_entrenamiento, ]

Para elaborar el modelo de distribución se utiliza el paquete dismo, el cual implementa varios algoritmos de modelado de distribución de especies. Cada algoritmo recibe como entradas un conjunto de datos de presencia de una especie y un conjunto de variables climáticas (en formato raster).

11.5.2 Modelo Maxent

El algoritmo Maxent se basa en el principio de máxima entropía, que sugiere que, con información incompleta, la mejor opción es la distribución menos sesgada posible que aún sea consistente con los datos. En este caso, Maxent encuentra la distribución de la especie que maximiza la entropía (es decir, la más uniforme posible) dado un conjunto de restricciones.

La implementación de Maxent que se incluye en el paquete dismo fue programada en el lenguaje Java. Por lo tanto, para ejecutarla es necesario instalar:

El Java Development Kit (JDK), para ejecutar aplicaciones desarrolladas en Java.
El paquete rJava para acceder a aplicaciones Java desde R.

11.5.2.1 Generación

# Los datos de clima deben convertirse al formato que usa el paquete raster
# debido a es este el que acepta el paquete dismo
clima <- raster::stack(clima)

# Ejecutar el modelo
modelo_maxent <- maxent(x = clima, p = entrenamiento)

# Aplicar el modelo entrenado a las variables climáticas 
# para generar un mapa de idoneidad del hábitat
prediccion <- predict(modelo_maxent, clima)

11.5.2.2 Evaluación

El siguiente bloqu de código genera métricas de desempeño del modelo utilizando los valores predichos en puntos de presencia y ausencia.

# terra::extract() extrae los valores del raster de predicción 
# en las coordenadas de evaluación
# eval_pres almacena los valores de idoneidad predichos 
# en los puntos de evaluación de presencia
eval_pres <- terra::extract(
  prediccion, 
  evaluacion[, c('decimalLongitude', 'decimalLatitude')]
)

# Generar puntos aleatorios dentro del área de estudio definida. 
# Estos puntos se asumen como ausencias de la especie.
ausencias <- randomPoints(mask = clima, n = 1000)

# eval_aus almacena los valores de idoneidad predichos
# en los puntos de ausencia
eval_aus <- terra::extract(
  prediccion, 
  ausencias
)

# Generar estadísticas de evaluación del modelo
resultado_evaluacion <- evaluate(p = eval_pres, a = eval_aus)

Las métricas de desempeño generadas en el bloque de código anterior se basan en la curva ROC y el AUC. La curva ROC (Receiver Operating Characteristic) es una herramienta gráfica que se utiliza para evaluar el rendimiento de un modelo de clasificación binaria, en este caso, para distinguir entre presencia y ausencia de una especie. El AUC (Area Under the Curve) es una métrica derivada de la curva ROC que resume su rendimiento en un solo valor.

Para interpretar la curva ROC y el AUC es necesario comprender los siguientes conceptos:

Verdaderos Positivos (True Positives, TP): son los casos en los que el modelo predice correctamente la presencia de la especie en un área donde la especie realmente está presente.
Falsos Positivos (False Positives, FP): son los casos en los que el modelo predice la presencia de la especie en un área donde la especie realmente no está presente.
Verdaderos Negativos (True Negatives, TN): son los casos en los que el modelo predice correctamente la ausencia de la especie en un área donde la especie realmente no está presente.
Falsos Negativos (False Negatives, FN): son los casos en los que el modelo predice la ausencia de la especie en un área donde la especie realmente está presente.

Las definiciones anteriores se resumen en la siguiente tabla:

	Predicción: Presencia	Predicción: Ausencia
Realidad: Presencia	Verdadero Positivo (TP)	Falso Negativo (FN)
Realidad: Ausencia	Falso Positivo (FP)	Verdadero Negativo (TN)

La tabla anterior es una matriz de confusión, una herramienta para evaluar el rendimiento de un modelo de clasificación. Resume los resultados de las predicciones del modelo comparándolas con los valores reales conocidos. En el contexto de los modelos de nicho ecológico, nos ayuda a entender cómo el modelo clasifica correctamente o incorrectamente las presencias y ausencias de una especie.

La curva ROC representa la relación entre la Tasa de Verdaderos Positivos (True Positive Rate, TPR) y la Tasa de Falsos Positivos (False Positive Rate, FPR) para diferentes umbrales de decisión del modelo.

TPR (Sensibilidad o recall): es la proporción de verdaderos positivos correctamente identificados por el modelo respecto al total de verdaderos positivos reales.

\[\text{TPR} = \frac{\text{Verdaderos Positivos}}{\text{Verdaderos Positivos} + \text{Falsos Negativos}}\]

FPR: es la proporción de falsos positivos identificados por el modelo respecto al total de negativos reales.

\[\text{FPR} = \frac{\text{Falsos Positivos}}{\text{Falsos Positivos} + \text{Verdaderos Negativos}}\]

Un umbral es un valor específico que se utiliza para convertir las predicciones continuas del modelo en clasificaciones binarias, es decir, en predicciones de presencia o ausencia de una especie. Los modelos como Maxent generan valores continuos que representan la idoneidad del hábitat o la probabilidad de presencia de una especie en cada punto en el espacio. Estos valores oscilan entre 0 (baja idoneidad) y 1 (alta idoneidad).

Para muchas aplicaciones prácticas, es necesario convertir estos valores continuos en categorías binarias de presencia o ausencia. Es necesario entonces definir un umbral (i.e. “un valor de corte”) que determina a partir de qué punto consideramos que la especie está presente.

Si el valor predicho es mayor o igual que el umbral, se predice presencia.
Si el valor predicho es menor que el umbral, se predice ausencia.

Es importante tener el cuenta el efecto que el umbral tiene en las predicciones. Si el umbral es alto, menos áreas serán clasificadas como presencia. Aumenta la especificidad (menos falsos positivos), pero puede disminuir la sensibilidad (más falsos negativos). Por otra parte, si el umbral es bajo, más áreas serán clasificadas como presencia. Aumenta la sensibilidad (menos falsos negativos), pero puede disminuir la especificidad (más falsos positivos).

Cada punto en la curva ROC representa un par (FPR, TPR) para un umbral específico. Al mover el umbral de decisión, cambiamos la sensibilidad y la especificidad del modelo, lo que se refleja en la curva. Una curva que se acerca más al punto superior izquierdo (TPR = 1, FPR = 0) indica un mejor rendimiento, ya que significa alta sensibilidad y baja tasa de falsos positivos. Una línea diagonal desde (0,0) hasta (1,1) representa un modelo que clasifica al azar (AUC = 0.5).

El AUC es el área bajo la curva ROC y proporciona una medida agregada del rendimiento del modelo en todos los umbrales posibles.

AUC = 1.0: rendimiento perfecto; el modelo clasifica correctamente todas las instancias positivas y negativas.
AUC = 0.5: rendimiento aleatorio; el modelo no tiene capacidad discriminativa.
AUC < 0.5: rendimiento peor que el azar; indica que el modelo está invertido en su capacidad predictiva.

El siguiente bloque de código muestra la curva ROC y el AUC para el modelo que se generó.

Código

# Datos para graficar la curva ROC
datos_roc <- data.frame(
  FPR = resultado_evaluacion@FPR,
  TPR = resultado_evaluacion@TPR,
  Umbral = resultado_evaluacion@t
)

# Valor AUC
auc <- resultado_evaluacion@auc

# Gráfico ggplot2
grafico_ggplot2 <-
  ggplot(
    datos_roc, 
    aes(
      x = FPR, 
      y = TPR,
      u = Umbral
    )
  ) +
  geom_line(
    color = "blue", 
    size = 1
  ) +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "red") +
  labs(title = paste("Curva ROC (AUC =", round(auc, 3), ")"),
       x = "Tasa de falsos positivos (FPR)",
       y = "Tasa de verdaderos positivos (TPR)") +
  theme_minimal()

# Gráfico plotly
ggplotly(grafico_ggplot2) |> 
  config(locale = 'es')

El valor de AUC es 0.953, lo que indica un excelente rendimiento del modelo. La curva ROC se eleva rápidamente cerca del eje y, lo que indica que el modelo tiene una alta tasa de verdaderos positivos incluso con bajos niveles de falsos positivos. Este comportamiento es característico de un buen modelo de clasificación.

11.5.2.3 Mapas

Seguidamente se presenta el mapa de idoneidad del hábitat.

Código

# Paleta de colores del modelo
colores_modelo <- colorNumeric(
  palette = c("white", "black"),
  values(prediccion),
  na.color = "transparent"
)

# Mapa
leaflet() |>
  addTiles(group = "Mapa general") |>
  addProviderTiles(
    providers$Esri.WorldImagery, 
    group = "Imágenes satelitales"
  ) |>  
  addProviderTiles(
    providers$CartoDB.Positron, 
    group = "Mapa blanco"
  ) |>
  addRasterImage( # capa raster de temperatura
    clima$wc2.1_10m_bio_1,
    colors = colores_temperatura, # paleta de colores
    opacity = 0.6,
    group = "Temperatura",
  ) |>
  addRasterImage( # capa raster de precipitación
    clima$wc2.1_10m_bio_12,
    colors = colores_precipitacion, # paleta de colores
    opacity = 0.6,
    group = "Precipitación",
  ) |>
  addRasterImage( # capa raster del modelo de distribución
    prediccion,
    colors = colores_modelo,
    opacity = 0.6,
    group = "Modelo de distribución",
  ) |>  
  addCircleMarkers(
    # capa de registros de presencia (puntos)
    data = presencia,
    stroke = F,
    radius = 3,
    fillColor = 'red',
    fillOpacity = 1,
    popup = paste(
      paste0("<strong>País: </strong>", presencia$country),
      paste0("<strong>Localidad: </strong>", presencia$locality),
      paste0("<strong>Fecha: </strong>", presencia$eventDate),
      paste0("<strong>Fuente: </strong>", presencia$institutionCode),
      paste0("<a href='", presencia$occurrenceID, "'>Más información</a>"),
      sep = '<br/>'
    ),
    group = "Registros de Bradypus variegatus"
  ) |>  
  addLegend(
    title = "Temperatura",
    values = values(clima$wc2.1_10m_bio_1),
    pal = colores_temperatura,
    position = "bottomleft",
    group = "Temperatura"
  ) |>
  addLegend(
    title = "Precipitación",
    values = values(clima$wc2.1_10m_bio_12),
    pal = colores_precipitacion,
    position = "bottomleft",
    group = "Precipitación"
  ) |>
  addLegend(
    title = "Modelo de distribución",
    values = values(prediccion),
    pal = colores_modelo,
    position = "bottomright",
    group = "Modelo de distribución"
  ) |>  
  addLayersControl(
    # control de capas
    baseGroups = c("Mapa general", "Imágenes satelitales", "Mapa blanco"),
    overlayGroups = c(
      "Temperatura",
      "Precipitación",
      "Modelo de distribución",
      "Registros de Bradypus variegatus"
    )
  ) |>
  hideGroup("Temperatura") |>
  hideGroup("Precipitación")

El siguiente mapa muestra las áreas en las que se predice presencia y ausencia de acuerdo con un umbral.

Código

# Definir el umbral
umbral <- 0.5

# Crear el raster binario
prediccion_binaria <- (prediccion >= umbral) * 1

# Crear la paleta de colores para el raster binario
colores_prediccion_binaria <- colorFactor(
  palette = c("transparent", "blue"),  # "transparent" para las áreas no adecuadas
  domain = c(0, 1),
  na.color = "transparent"
)

# Mapa
leaflet() |>
  addTiles(group = "Mapa general") |>
  addProviderTiles(
    providers$Esri.WorldImagery, 
    group = "Imágenes satelitales"
  ) |>
  addProviderTiles(
    providers$CartoDB.Positron, 
    group = "Mapa blanco"
  ) |>
  addRasterImage(
    prediccion_binaria,
    colors = colores_prediccion_binaria,
    opacity = 0.6,
    group = "Modelo de distribución binario",
  ) |>
  addCircleMarkers(
    data = presencia,
    stroke = FALSE,
    radius = 3,
    fillColor = 'red',
    fillOpacity = 1,
    popup = paste(
      paste0("<strong>País: </strong>", presencia$country),
      paste0("<strong>Localidad: </strong>", presencia$locality),
      paste0("<strong>Fecha: </strong>", presencia$eventDate),
      paste0("<strong>Fuente: </strong>", presencia$institutionCode),
      paste0("<a href='", presencia$occurrenceID, "'>Más información</a>"),
      sep = '<br/>'
    ),
    group = "Registros de Bradypus variegatus"
  ) |>
  addLegend(
    title = "Modelo de distribución binario",
    labels = c("Ausencia", "Presencia"),
    colors = c("transparent", "blue"),
    position = "bottomright",
    group = "Modelo de distribución binario"
  ) |>
  addLayersControl(
    baseGroups = c("Mapa general", "Imágenes satelitales", "Mapa blanco"),
    overlayGroups = c(
      "Modelo de distribución binario",
      "Registros de Bradypus variegatus"
    )
  )

11.6 Recursos de interés

Phillips, S. J., Anderson, R. P., & Schapire, R. E. (2006). Maximum entropy modeling of species geographic distributions. Ecological Modelling, 190(3-4), 231-259. https://doi.org/10.1016/j.ecolmodel.2005.03.026