Modelos de distribución de especies

Práctica 8, Datos Espaciales en Biodiversidad

Máster en Áreas Protegidas, Recursos Naturales y Biodiversidad

Profesores

María Victoria Jiménez Franco

José Francisco Calvo Sendín

Departamento

Ecología e Hidrología, Facultad de Biología

Fecha de publicación

9 de abril de 2024

Fecha de modificación

14 de abril de 2024

Introducción

En esta práctica se realiza una introducción a los modelos de distribución de especies (species distribution models, SDM), centrada en diferentes aspectos del procedimiento, como la preparación de los datos, el uso de diversos algoritmos, la proyección y evaluación de los modelos resultantes, y el ensamble de diferentes modelos. Realizaremos ejercicios para modelar datos de presencia/ausencia y datos únicamente de presencias (con generación de pseudausencias). Para cada uno de los casos, en primer lugar se prepararán los datos mediante procedimientos de manipulación de datos vectoriales y ráster similares a los utilizados en prácticas anteriores de la asignatura.

Preparación

Los datos, mapas y funciones necesarios para el desarrollo de la práctica se encuentran disponibles en el archivo RData de la asignatura. Una vez iniciado R, cargaremos este archivo desde el servidor:

load( url( "http://webs.um.es/jfcalvo/deb.RData" ) )

Usando ls() podemos ver los objetos cargados. Utiliza la función info para obtener información sobre ellos; por ejemplo: info( "taxiados" ) o info( "spain" ). Teclea info() para más información.

Para realizar esta práctica se necesita la instalación adicional de los paquetes de R terra, geodata, biomod2 y rgbif. Si no los tenemos instalados usaremos:

install.packages( c( "terra", "geodata", "biomod2", "rgbif" ) )

Una vez instalados los cargaremos en memoria:

library( terra )
library( geodata )
library( biomod2 )
library( rgbif )

Descomprimiremos también los mapas que necesitaremos en la práctica, usando la función unwrap y cambiaremos el sistema de referencia de coordenadas del mapa murcia:

unwrap( murcia ) -> murcia
unwrap( taxiados ) -> taxiados
unwrap( mdt200 ) -> mdt200
unwrap( corine90 ) -> corine90
unwrap( spain ) -> spain
project( murcia, crs( mdt200 ) ) -> murcia

Modelos de distribución con presencias y ausencias

Preparación de los datos

Trabajaremos con los datos del objeto taxiados en los que está registrada la presencia/ausencia de varias especies de aves en 140 itinerarios de censo realizados en sierras de la Región de Murcia. Para este ejercicio utilizaremos los datos del mito (Aegithalos caudatus).

plot( murcia )
plot( taxiados, "Aegithalos.caudatus", add = TRUE )

Inicialmente necesitamos extraer los datos (presencias y ausencias) del mito, y las coordenadas de cada taxiado. Como taxiados es un objeto de líneas consideraremos los centroides:

aegcau <- taxiados$Aegithalos.caudatus
xy <- geom( centroids( taxiados ) )[ , 3 : 4 ]

Como variables ambientales utilizaremos un mapa ráster con capas de altitud, temperatura media anual, precipitación anual y porcentaje de superficie de bosque. Para crearlo seguiremos un procedimiento similar al utilizado en la práctica 4 de la asignatura. En primer lugar, descargaremos los mapas bioclimáticos con la función worldclim_country del paquete geodata:

worldclim_country( country = "Spain" , var = "bio", path = "./descargas/" ) -> spain_bio
project( spain_bio, crs( mdt200 ) ) -> spain_bio
resample( spain_bio, mdt200 ) -> murcia_bio

A continuación rasterizaremos el conjunto de los tipos de uso “bosque” del mapa CORINE, estimando su proporción en cada cuadrícula (argumento cover = TRUE y opción background = 0):

subset( corine90, corine90$CODE_90 > 300 & corine90$CODE_90 < 320 ) -> bosque
rasterize( bosque, mdt200, cover = TRUE, background = 0 ) -> rbosque

Finalmente crearemos el mapa ráster, añadiendo también el mapa de elevaciones, y lo recortaremos con los límites territoriales de la Región de Murcia:

amb <- c( mdt200, murcia_bio[[ 1 ]], murcia_bio[[ 12 ]], rbosque )
crop( amb, murcia, mask = TRUE, snap = "out" ) -> amb
names( amb ) <- c( "alt", "temp", "prec", "bosq" )
plot( amb )

Ahora ya estamos en disposición de formatear los datos biológicos y ambientales para realizar los SDM. Utilizaremos la función BIOMOD_FormatingData, proporcionando los datos de presencia/ausencia de la especie, las coordenadas de las unidades de muestreo y el mapa de variables ambientales:

sdmData <- BIOMOD_FormatingData( resp.var = aegcau, resp.xy =  xy, expl.var = amb, resp.name = "Mito" )
sdmData

Modelado

Para realizar los modelos utilizaremos la función BIOMOD_Modeling. Podemos ver los diferentes algoritmos de modelado disponibles, además de las múltiples opciones de la función, con:

?BIOMOD_Modeling

Elegiremos inicialmente, por ejemplo, el modelo GLM (Generalized Linear Model). También utilizaremos inicialmente la totalidad de los datos, prescindiendo así de realizar un procedimiento de validación. Esta opción puede modificarse con el argumento CV.perc, que permite especificar la proporción de datos a utilizar en los modelos. La función que realiza el modelo es BIOMOD_Modeling:

sdmGLM <- BIOMOD_Modeling( bm.format = sdmData, models = "GLM", CV.perc = 1 )

La proyección de los modelos se realiza en función de los datos ambientales, que podrían extenderse a áreas espaciales o escalas temporales diferentes. En nuestro caso, realizaremos la proyección sobre el área de estudio bajo las condiciones actuales (es decir, utilizando las variables ambientales del objeto amb). La función que realiza este proceso es BIOMOD_Projection:

sdmGLMProj <- BIOMOD_Projection( bm.mod = sdmGLM, new.env = amb, proj.name = "Mito")

La visualización del mapa resultante se obtiene con la función plot aplicada sobre la función get_predictions.

plot( get_predictions( sdmGLMProj ) )

Con la función get_evaluations se obtienen diferentes métricas de evaluación:

get_evaluations( sdmGLM )

Finalmente, la función bm_PlotResponseCurves proporciona las curvas de respuesta de las variables ambientales:

bm_PlotResponseCurves( sdmGLM )

Ejercicio 1

Realizaremos en este ejercicio un SDM con los datos de presencia/ausencia del mito y el método Generalized Additive Model):

sdmGAM <- BIOMOD_Modeling( bm.format = sdmData, models = "GAM", CV.perc = 1 )
sdmGAMProj <- BIOMOD_Projection( bm.mod = sdmGAM, new.env = amb, proj.name = "Mito")
plot( get_predictions( sdmGAMProj ) )

get_evaluations( sdmGAM )
bm_PlotResponseCurves( sdmGAM )

Ensambles

La función BIOMOD_Modeling permite realizar simultáneamente varios modelos (diferentes algoritmos), con la posibilidad de realizar un ensamble, o modelado conjunto final con todos o una selección de ellos. Siguiendo con el ejemplo del mito, elaboraremos diversos modelos, considerando en este caso una sola réplica de cada uno (con todos los datos) para simplificar la visualización de los resultados. Dado que el ensamble se realiza en función de las métricas de evaluación, seleccionaremos únicamente ROC:

sdm <- BIOMOD_Modeling( bm.format = sdmData, models = c( "GLM", "GBM", "GAM", "ANN", "FDA", "MARS", "RF" ), CV.perc = 1, metric.eval = "ROC" )
sdmProj <- BIOMOD_Projection( bm.mod = sdm, new.env = amb, proj.name = "Mito")
plot( get_predictions( sdmProj ) )

Si quisieramos representar el resultado de un único modelo utilizaríamos, por ejemplo:

plot( get_predictions( sdmProj )[[ 7 ]] )

La evaluación de los modelos y las curvas de respuesta las obtenemos con:

get_evaluations( sdm )
bm_PlotResponseCurves( sdm )

El ensamble lo haremos con la función BIOMOD_EnsembleModeling y la proyección con BIOMOD_EnsembleForecasting:

sdmE <- BIOMOD_EnsembleModeling( bm.mod = sdm ) 
sdmEProj <- BIOMOD_EnsembleForecasting( bm.em = sdmE, new.env = amb, proj.name = "Mito" )
plot( get_predictions( sdmEProj ) )

La evaluación del modelo final y las curvas de respuesta las obtenemos con:

get_evaluations( sdmE )
bm_PlotResponseCurves( sdmE )

Validación

Los procedimientos de modelado SDM incluyen habitualmente la aplicación de métodos validación consistentes en seleccionar (habitualmente de forma aleatoria) una proporción de los datos para el calibrado del modelo y el resto para su validación. En estos casos, por tanto, se requiere la realización de varias réplicas. Necesitaremos utilizar los argumentos CV.nb.rep (para especificar el número de réplicas a realizar) y CV.perc (para especificar la proporción de datos a utilizar en los modelos). Además utilizaremos otro argumento (var.import) que especifica un número de permutaciones que utilizará la función para obtener un índice de importancia de cada una de las variables ambientales utilizadas en el modelo. Analizaremos este procedimiento de nuevo con los datos de presencia/ausendia del mito, y seleccionando en este caso, por ejemplo, el algoritmo Random Forest:

sdmRF <- BIOMOD_Modeling( bm.format = sdmData, models = "RF", CV.nb.rep = 6, CV.perc = 0.7, var.import = 10 )
sdmRFProj <- BIOMOD_Projection( bm.mod = sdmRF, new.env = amb, proj.name = "Mito")
plot( get_predictions( sdmRFProj ) )

Como tenemos varias réplicas del modelo, para representar el “modelo promedio” podríamos utilizar:

plot( mean( get_predictions( sdmRFProj ) ) )

En este caso, dado que tenemos varías réplicas con sus corrrespondientes métricas de evaluación, resulta apropiado usar la función bm_PlotEvalBoxplot:

bm_PlotEvalBoxplot( sdmRF, group.by = c( "algo", "algo" ) )

Por otra parte, la función bm_PlotVarImpBoxplot proporciona la información sobre la importancia de las variables ambientales:

bm_PlotVarImpBoxplot( sdmRF, group.by = c( "expl.var", "algo", "algo" ) )

Finalmente, la función bm_PlotResponseCurves proporciona las curvas de respuesta de las variables ambientales:

bm_PlotResponseCurves( sdmRF )

Modelos de distribución con presencias y pseudoausencias

Los modelos de distribución de especies se pueden realizar también con información únicamente de presencias, que suelen corresponder a registros o citas de especies procedentes, por ejemplo, de programas de ciencia ciudadana. El objetivo de esta sección es realizar el modelado de datos de presencias para elaborar mapas de idoneidad de hábitat en función de diferentes variables ambientales.

Realizaremos un ejercicio con datos de la alondra ricotí (Chersophilus duponti) obtenidos a partir de la base de datos mundial del Sistema Global de Información sobre Biodiversidad (GBIF, https://www.gbif.org/). Utilizaremos datos de localización de esta especie en España entre los años 2020 y 2024 (que obtendremos siguiendo un procedimiento similar al realizado en la práctica 2 de la asignatura):

occ_data( scientificName = "Chersophilus duponti", country = "ES", year = "2020,2024", limit = 1000 ) -> chedup_ES
vect( chedup_ES$data, geom = c( "decimalLongitude", "decimalLatitude" ), crs("LonLat") ) -> chedup
project( chedup, crs( spain ) ) -> chedup
chedup
plot( spain )
plot( chedup, add = TRUE, cex = 0.5, col = 4 )

A continuación generaremos un vector de presencias y las coordenadas de las localizaciones:

chedup$pres <- 1
xy2 <- geom( chedup )[ , 3 : 4 ]

Como variables ambientales utilizaremos las variables bioclimáticas ya descargadas (objeto spain_bio), que recortaremos con los límites territoriales de España peninsular y Baleares:

spain_bio <- crop( spain_bio, spain, mask = TRUE )
plot( spain_bio )

Por último, formateamos los datos para aplicar las funciones de biomod2. En este caso hay que utilizar los argumentos PA.nb.rep y PA.strategy para generar las pseudoausencias necesarias para la realización de los modelos. También es importante utilizar el argumento filter.raster = TRUE para eliminar duplicados (presencias coincidentes en una misma celda del mapa ráster de variables ambientales):

sdmData2 <- BIOMOD_FormatingData( resp.var = chedup$pres, resp.xy =  xy2, expl.var = spain_bio, resp.name = "Alondra.ricotí", PA.nb.rep = 1, PA.strategy = "random", filter.raster = TRUE )
sdmData2

Como ejemplo de modelo que utiliza pseudoausencias utilizaremos el algoritmo MAXNET, que realiza un modelo de Máxima Entropía:

sdmME <- BIOMOD_Modeling( bm.format = sdmData2, models = "MAXNET", CV.nb.rep = 6, CV.perc = 0.7, var.import = 10 )
sdmMEProj <- BIOMOD_Projection( bm.mod = sdmME, new.env = spain_bio, proj.name = "Alonda.ricotí" )
plot( mean( get_predictions( sdmMEProj ) ) )

La evaluación del modelo la obtenemos con:

bm_PlotEvalBoxplot( sdmME, group.by = c( "algo", "algo" ) )

La información sobre la importancia de las variables bioclimáticas la obtenemos con:

bm_PlotVarImpBoxplot( sdmME, group.by = c( "expl.var", "algo", "algo" ) )

Las curvas de respuesta las obtenemos con:

bm_PlotResponseCurves( sdmME )

Evaluación

Realiza la prueba de evaluación disponible en la herramienta Exámenes del Aula Virtual.

Bibliografía

Guisan A, Thuiller W., & Zimmermann, N. E. (2017) Habitat suitability and distribution models: with applications in R. Cambridge University Press.
Hijmans R (2024) _terra: Spatial Data Analysis_. R package version 1.7-71, https://CRAN.R-project.org/package=terra
Thuiller W, Georges D, Gueguen M, Engler R, Breiner F, Lafourcade B, Patin R

_biomod2: Ensemble Platform for Species Distribution Modeling_. R package version 4.2-4, https://CRAN.R-project.org/package=biomod2.

Descripción de los mapas y funciones

Utiliza la función info. Por ejemplo: info( "taxiados" ) o info( "spain" ).