Transferencias de modelos de nicho ecológico
KUENM, MOP, procedimientos de transferencia
Kevin López-Reyes
Universidad Nacional Autonoma de México
¿Qué es una transferencia?
En MNE, transferir un modelo implica proyectar las relaciones especie–ambiente estimadas en una region y periodo especifico hacia nuevos escenarios.
Implica dos procesos en el espacio ambiental:
Interpolación : transferencias hacia condiciones ambientales dentro de las utilizadas para calibrar el modelo.Extrapolación : transferencias hacia condiciones ambientales fuera de las utilizadas para calibrar el modelo.
Tipos de transferencia
Las transferencias pueden ocurrir a traves de:
el espacio : una región geográfica distinta
el tiempo : pasado o futuro;
ambos : otra región en otro tiempo
¿Por qué quisieran transferir un modelo de nicho?
Invasiones biológicas
Sapo de caña (Rhinella marina )
El modelo se calibró con las condiciones climáticas del rango nativo de la especie en América y se transfirió a Australia, donde fue introducida.
Origen: Centroamérica y Sudamérica.Destino: Australia.Pregunta: ¿Qué regiones de Australia presentan condiciones climáticas adecuadas?
¿Por qué quisieran transferir un modelo de nicho?
Proyecciones a futuro
Panda gigante (Ailuropoda melanoleuca )
Un modelo Maxent calibrado con registros y condiciones ambientales actuales fue transferido a escenarios de 2050 y 2070.
Origen: condiciones ambientales actuales en Sichuan, China.Destino: condiciones proyectadas para 2050 y 2070.Pregunta: ¿Cómo cambiará la distribución del panda gigante bajo distintos escenarios climáticos?
Liu et al. (2022) . Bajo RCP8.5, el modelo proyectó una pérdida neta de hábitat adecuado de 25% para 2050 y 38% para 2070.
Objetivo de la práctica
Desarrollar un modelo de nicho ecológico con kuenm2 , desde la limpieza y preparación de los datos hasta la calibración, selección, y transferencia del modelo.
Durante las transferencias compararemos:
extrapolación libre ;clamping ;truncamiento o no extrapolación ;riesgo de extrapolación mediante el análisis MOP .
¿Qué es kuenm2?
kuenm2 es un paquete de R para desarrollar modelos de nicho ecológico de manera detallada y reproducible.
Actualmente permite trabajar con:
Maxnet , basado en máxima entropía;GLM , modelos lineales generalizados.
Incluye herramientas para preparar y explorar los datos, calibrar y seleccionar modelos, realizar transferencias y analizar su incertidumbre.
¿Qué podemos hacer con kuenm2?
Antes de modelar
limpiar y preparar registros de presencia;
explorar los datos geográfica y ambientalmente;
crear particiones de entrenamiento y evaluación.
Durante el modelado
calibrar múltiples modelos candidatos;
seleccionar modelos mediante varios criterios;
examinar importancia de variables y curvas de respuesta;
proyectar los modelos hacia escenarios nuevos.
Después de modelar
comparar predicciones;
cuantificar su variabilidad;
analizar el riesgo de extrapolación mediante MOP.
Flujo de trabajo de kuenm2
Estructura de la práctica
Seguiremos las etapas descritas en las guías oficiales:
Limpieza básica de los datos.
Preparación de los datos de calibración.
Calibración y selección de modelos.
Ajuste y exploración de los modelos seleccionados.
Predicción hacia un escenario individual.
Análisis de incertidumbre.
Fuente de esta práctica
Esta práctica es una adaptación de las guías oficiales de kuenm2, desarrolladas por Weverton C. F. Trindade, Luis F. Arias-Giraldo, Luis Osorio-Olvera, A. Townsend Peterson y Marlon E. Cobos.
La documentación original contiene más opciones. Recomiendo revisarla!
Documentación original
El código de kuenm2 se distribuye bajo licencia GPL ≥ 3 .
Sistema de estudio
Lagartija espinosa azul yucateca (Sceloporus serrifer ). Una especie de lagartija de la Peninsula de Yucatan.
Calibrar un modelo en el presente
Transferir el modelo al periodo 2040-2060
Revisar los resultados
Analizar el riesgo de extrapolación
Insumos
Descargar los datos en el siguiente enlace: DESCARGAR
Cargamos los paquetes necesarios para la practica
library (kuenm2)library (terra)library (leaflet)library (magrittr)
1. Limpieza básica de los datos
Antes de estimar un modelo debemos limpiar nuestros datos!
En esta sección:
inspeccionaremos los registros originales;
aplicaremos distintos criterios de limpieza;
compararemos los datos antes y después del proceso.
Cargar los registros de presencia
<- read.csv ("occ_joint.csv" )head (occ_raw)
id uuid observed_on_string
1 1276590 653d1567-91c5-48d4-82ea-2873069d6877 2012-05-04
2 2349231 56a61b7f-9c97-4416-bb2e-5df72d44449e 2011-06-19
3 2695920 3ce7ede8-ec7a-496b-8b54-f3dea926f666 2015-11-29 9:36:42 AM CST
4 3697300 490d01d7-5807-4dc9-8117-498141bea807 2007/03/16 4:12 PM HST
5 3697301 93c6d784-6764-402e-bb4e-2aacc70ce16b 2007/03/16 4:13 PM HST
6 3697302 665bc265-1b38-41c9-88cb-c08712127976 2007/03/16 4:13 PM HST
observed_on time_observed_at time_zone user_id
1 2012-05-04 Mexico City 1315
2 2011-06-19 Mexico City 77824
3 2015-11-29 2015-11-29 15:36:42 UTC Central Time (US & Canada) 16074
4 2007-03-16 2007-03-17 02:12:00 UTC Hawaii 238859
5 2007-03-16 2007-03-17 02:13:00 UTC Hawaii 238859
6 2007-03-16 2007-03-17 02:13:00 UTC Hawaii 238859
user_login user_name created_at
1 escalante-pasos Jorge Armín Escalante-Pasos 2015-03-07 17:07:29 UTC
2 adorantes Alfredo Dorantes Euan 2015-11-04 22:23:15 UTC
3 sergiofernando121 Sergio Fernando 2016-02-19 02:49:41 UTC
4 rafas-lara Rafa Lara 2016-07-18 07:54:17 UTC
5 rafas-lara Rafa Lara 2016-07-18 07:54:18 UTC
6 rafas-lara Rafa Lara 2016-07-18 07:54:18 UTC
updated_at quality_grade license
1 2021-11-19 01:07:31 UTC research CC-BY
2 2021-10-13 14:18:57 UTC research
3 2021-10-13 14:19:00 UTC research CC-BY-NC
4 2022-10-12 04:30:26 UTC research
5 2022-10-12 04:30:26 UTC research
6 2022-10-12 04:30:26 UTC research
url
1 http://conabio.inaturalist.org/observations/1276590
2 http://conabio.inaturalist.org/observations/2349231
3 http://www.inaturalist.org/observations/2695920
4 http://www.inaturalist.org/observations/3697300
5 http://www.inaturalist.org/observations/3697301
6 http://www.inaturalist.org/observations/3697302
image_url
1 https://inaturalist-open-data.s3.amazonaws.com/photos/1601390/medium.jpg
2 https://static.inaturalist.org/photos/2609829/medium.JPG
3 https://inaturalist-open-data.s3.amazonaws.com/photos/3034176/medium.jpeg
4 https://static.inaturalist.org/photos/4291334/medium.jpeg
5 https://static.inaturalist.org/photos/4291336/medium.jpeg
6 https://static.inaturalist.org/photos/4291337/medium.jpeg
sound_url tag_list description
1 NA Macho subadulto
2 NA Reptiles Yucatán, Reptiles Homún
3 NA
4 NA
5 NA
6 NA
num_identification_agreements num_identification_disagreements
1 3 0
2 2 0
3 2 0
4 1 0
5 1 0
6 1 0
captive_cultivated oauth_application_id
1 false NA
2 false NA
3 false 2
4 false NA
5 false NA
6 false NA
place_guess latitude longitude
1 San Antonio Kaua 20.94463 -89.58681
2 Homun, Yucatán 20.73728 -89.28595
3 Calle 85 30, Santa Isabel, 97370 Kanasín, Yuc., México 20.91694 -89.57487
4 Ocosingo, Chis., México 16.89572 -90.96684
5 Ocosingo, Chis., México 16.89778 -90.96877
6 Ocosingo, Chis., México 16.89687 -90.96877
positional_accuracy private_place_guess private_latitude private_longitude
1 31 NA NA NA
2 15 NA NA NA
3 33 NA NA NA
4 61 NA NA NA
5 61 NA NA NA
6 61 NA NA NA
public_positional_accuracy geoprivacy taxon_geoprivacy coordinates_obscured
1 31 open false
2 15 open false
3 33 open false
4 61 open false
5 61 open false
6 61 open false
positioning_method positioning_device species_guess
1 Sceloporus serrifer serrifer
2 Lagartija espinosa azul
3 gps gps Yucatecan Blue Spiny Lizard
4 Lagartija espinosa azul
5 Lagartija espinosa azul
6 Lagartija espinosa azul
scientific_name common_name iconic_taxon_name taxon_id
1 Sceloporus serrifer Yucatecan Blue Spiny Lizard Reptilia 36178
2 Sceloporus serrifer Yucatecan Blue Spiny Lizard Reptilia 36178
3 Sceloporus serrifer Yucatecan Blue Spiny Lizard Reptilia 36178
4 Sceloporus serrifer Yucatecan Blue Spiny Lizard Reptilia 36178
5 Sceloporus serrifer Yucatecan Blue Spiny Lizard Reptilia 36178
6 Sceloporus serrifer Yucatecan Blue Spiny Lizard Reptilia 36178
Seleccionar y renombrar columnas
<- occ_raw[, c ("scientific_name" ,"longitude" ,"latitude" )]names (occ) <- c ("species" , "longitude" , "latitude" )head (occ)
species longitude latitude
1 Sceloporus serrifer -89.58681 20.94463
2 Sceloporus serrifer -89.28595 20.73728
3 Sceloporus serrifer -89.57487 20.91694
4 Sceloporus serrifer -90.96684 16.89572
5 Sceloporus serrifer -90.96877 16.89778
6 Sceloporus serrifer -90.96877 16.89687
Cargar variables ambientales
<- rast (list.files ("M_var/set_1" ,pattern = ".tif" ,full.names = TRUE ))plot (variables)
Definir el área de calibración
Construimos un buffer de 200 km alrededor de los registros y recortamos las variables ambientales.
# Convertir los registros en un objeto espacial <- vect (occ, geom = c ("longitude" , "latitude" ),crs = "EPSG:4326" )# Proyectar a coordenadas en metros <- project (registros, "EPSG:32616" )# Crear un buffer de 200 km <- buffer (registros_utm, width = 200000 )# Unir los buffers <- aggregate (M_buffer)# Regresar al CRS de las variables <- project (M_buffer, crs (variables))# Recortar las variables <- crop (variables, M_buffer)<- mask (M_var, M_buffer)plot (M_var)
¿Qué estamos asumiendo al definir M?
Los ambientes que han sido accesibles para la especie, dónde se generan los puntos de fondo
Si aumentamos el buffer a 500 km, ¿cambiará el modelo aunque las presencias sean las mismas?
Sí. Cambiarían los ambientes de fondo y, por lo tanto, el contraste utilizado para estimar las relaciones especie–ambiente.
Limpieza básica de los registros
<- initial_cleaning (data = occ,species = "species" ,x = "longitude" ,y = "latitude" ,remove_na = TRUE ,remove_empty = TRUE ,remove_duplicates = TRUE ,by_decimal_precision = TRUE ,decimal_precision = 2 )nrow (occ)
Eliminar registros duplicados por celda
Conservamos una sola presencia dentro de cada celda de las variables ambientales.
set.seed (123 )<- remove_cell_duplicates (data = occ_clean,x = "longitude" ,y = "latitude" ,raster_layer = M_var[[1 ]])nrow (occ_clean)
Registros antes y después de la limpieza
2. Preparación de los datos
prepare_data():
extrae los valores ambientales de las presencias;
genera puntos de fondo dentro del área de calibración;
divide los registros para evaluar los modelos;
crea las combinaciones de parámetros que serán calibradas.
2.1. Configuración
Utilizaremos:
algoritmo: Maxnet ;
cuatro particiones;
1,000 puntos de fondo;
dos tipos de respuesta;
tres valores de regularización.
Esto produce seis modelos candidatos .
set.seed (123 )<- prepare_data (algorithm = "maxnet" ,occ = occ_clean_cell,x = "longitude" ,y = "latitude" ,raster_variables = M_var,min_number = 4 ,species = "Sceloporus serrifer" ,partition_method = "kfolds" ,n_partitions = 4 ,n_background = 1000 ,features = c ("l" , "q" ),r_multiplier = c (0.5 , 1 , 2 ))
2.2. Datos preparados
prepared_data object summary
============================
Species: Sceloporus serrifer
Number of Records: 1053
- Presence: 53
- Background: 1000
Partition Method: kfolds
- Number of kfolds: 4
Continuous Variables:
- bio_1, bio_12, bio_15, bio_4
Categorical Variables: None
PCA Information: PCA not performed
Weights: No weights provided
Calibration Parameters:
- Algorithm: maxnet
- Number of candidate models: 6
- Features classes (responses): l, q
- Regularization multipliers: 0.5, 1, 2
El resumen permite comprobar:
número de presencias y puntos de fondo;
variables utilizadas;
método de partición;
número de modelos candidatos.
3. Calibración de modelos
Los seis modelos candidatos serán evaluados utilizando cuatro criterios:
ausencia de curvas cóncavas;
tasa de omisión 10%;
significancia mediante ROC parcial;
menor complejidad mediante AICc.
Respuestas concavas
“Concave (or bimodal) curves indicate that higher suitability is found at variable values around a point of lower suitability”
No tienen tanto sentido biológico
3.1. Calibrar los modelos candidatos
<- calibration (data = datos_calibracion,remove_concave = TRUE ,error_considered = 10 ,omission_rate = 10 ,parallel = FALSE ,ncores = 1 )
|
| | 0%
|
|======================================================================| 100%
|
| | 0%
|
|======================= | 33%
|
|=============================================== | 67%
|
|======================================================================| 100%
|
| | 0%
|
|======================= | 33%
|
|=============================================== | 67%
|
|======================================================================| 100%
3.2. Modelos seleccionados
calibration_results object summary (maxnet)
=============================================================
Species: Sceloporus serrifer
Number of candidate models: 6
- Models removed because they failed to fit: 0
- Models identified with concave curves: 3
- Model with concave curves removed
- Models removed with non-significant values of pROC: 0
- Models removed with omission error > 10%: 0
- Models removed with delta AIC > 2: 0
Selected models: 3
- Up to 5 printed here:
ID Formulas Features R_multiplier
1 1 ~bio_1 + bio_12 + bio_15 + bio_4 -1 l 0.5
2 2 ~bio_1 + bio_12 + bio_15 + bio_4 -1 l 1.0
3 3 ~bio_1 + bio_12 + bio_15 + bio_4 -1 l 2.0
pval_pROC_at_10.mean Omission_rate_at_10.mean dAIC Parameters
1 0 0.0948 0.0000000 4
2 0 0.0948 0.0907178 4
3 0 0.0755 0.4388027 4
$ selected_models
ID Formulas R_multiplier Features
1 1 ~bio_1 + bio_12 + bio_15 + bio_4 -1 0.5 l
2 2 ~bio_1 + bio_12 + bio_15 + bio_4 -1 1.0 l
3 3 ~bio_1 + bio_12 + bio_15 + bio_4 -1 2.0 l
Omission_rate_at_10.mean Omission_rate_at_10.sd Mean_AUC_ratio_at_10.mean
1 0.0948 0.0395 1.206010
2 0.0948 0.0395 1.204976
3 0.0755 0.0027 1.201239
Mean_AUC_ratio_at_10.sd pval_pROC_at_10.mean pval_pROC_at_10.sd AICc
1 0.008649086 0 0 713.0060
2 0.006980445 0 0 713.0967
3 0.006494816 0 0 713.4448
Parameters Is_concave dAIC
1 4 FALSE 0.0000000
2 4 FALSE 0.0907178
3 4 FALSE 0.4388027
Los modelos seleccionados:
omiten menos del 10% de las presencias de evaluación;
predicen mejor que un modelo aleatorio;
presentan un delta AICc menor que 2.
4. Ajustar los modelos seleccionados
fit_selected() vuelve a ajustar las configuraciones seleccionadas utilizando:
cuatro réplicas correspondientes a las particiones;.
<- fit_selected (calibration_results = calibracion,replicate_method = "kfolds" ,n_replicates = 4 ,type = "cloglog" ,progress_bar = FALSE ,verbose = FALSE ,seed = 123 )
4.1. Modelos finales
fitted_models object summary
============================
Species: Sceloporus serrifer
Algortihm: maxnet
Number of fitted models: 3
Models fitted with 4 replicates
Cada configuración seleccionada contiene:
cuatro modelos (particiones);
un modelo completo ajustado con todos los datos.
4.2. Curvas de respuesta
Las curvas muestran cómo cambia la idoneidad estimada a lo largo de cada variable ambiental.
all_response_curves (models = modelos_finales,show_variability = TRUE ,show_lines = TRUE ,extrapolate = FALSE ,mfrow = c (2 , 2 ),col = "skyblue4" )
Las líneas verticales discontinuas delimitan los valores utilizados para ajustar los modelos.
4.3. Importancia de las variables
<- variable_importance (models = modelos_finales,progress_bar = FALSE ,verbose = FALSE )
plot_importance (main = "Importancia de las variables" )
5. Transferencia a un escenario futuro
Cargar las variables climáticas futuras
Utilizaremos las condiciones proyectadas por MIROC6 para el escenario SSP5-8.5 durante 2041–2060 .
<- rast ("G_var/set_1/future/bios.tif" )<- future_all[[c (1 , 4 , 12 , 15 )]]names (G_future) <- names (M_var)
class : SpatRaster
size : 1080, 2160, 4 (nrow, ncol, nlyr)
resolution : 0.1666667, 0.1666667 (x, y)
extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326)
source : bios.tif
names : bio_1, bio_12, bio_15, bio_4
min values : -52.8, 11.2, 0.0, 0.0
max values : 33.2, 2313.7, 8010.2, 236.3
Recortar las variables futuras
Las variables futuras se recortan utilizando la misma área empleada para calibrar el modelo.
<- crop (G_future, M_buffer)<- mask (G_var, M_buffer)<- mask (G_var, M_var[[1 ]])names (G_var) <- names (M_var)plot (G_var)
Verificar las variables
Los predictores presentes y futuros deben tener los mismos nombres
[1] "bio_1" "bio_12" "bio_15" "bio_4"
[1] "bio_1" "bio_12" "bio_15" "bio_4"
¿Por qué?
Recordemos que la transferencia es proyectar una relación especies-clima estimada hacia otro escenario (region o tiempo)
Esas relaciones están estimadas con unas variables especificas, por eso deben mantenerse iguales.
Proyección en las condiciones actuales
<- predict_selected (models = modelos_finales,new_variables = M_var,type = "cloglog" ,consensus = "mean" ,extrapolation_type = "NE" ,progress_bar = FALSE )<- prediccion_actual$ General_consensus$ meanplot (mapa_actual)
Idoneidad ambiental actual
Procedimientos de transferencia
Cuando las condiciones de la región de transferencia exceden los rangos ambientales de calibración, podemos utilizar tres procedimientos:
Extrapolación libre — E: Las respuestas estimadas continúan fuera del rango de calibración.
Extrapolación con clamping — EC: Los valores ambientales que exceden el rango se sustituyen por el límite mínimo o máximo conocido.
No extrapolación o truncamiento — NE: Las celdas con al menos una variable fuera del rango de calibración reciben idoneidad cero.
Dentro de los rangos de calibración, los tres procedimientos producen la misma predicción.
¿Por qué no simplemente incrementamos el area de calibración?
¿Por qué no simplemente incrementamos el area de calibración?
Si usamos una metrica dependiente de umbral, puede incrementar las metricas de evalución
Le estamos diciendo al modelo que la especie pudo muestrear esos ambientes y no está ahí porque no es favorable, lo cual no es cierto
Transferencia al futuro
Realizaremos la transferencia utilizando los tres procedimientos:
E: extrapolación libre.EC: extrapolación con clamping .NE: sin extrapolación.
<- predict_selected (models = modelos_finales,new_variables = G_var,type = "cloglog" ,consensus = "mean" ,extrapolation_type = "E" ,progress_bar = FALSE )<- predict_selected (models = modelos_finales,new_variables = G_var,type = "cloglog" ,consensus = "mean" ,extrapolation_type = "EC" ,progress_bar = FALSE )<- predict_selected (models = modelos_finales,new_variables = G_var,type = "cloglog" ,consensus = "mean" ,extrapolation_type = "NE" ,progress_bar = FALSE )
Obtener las predicciones
<- prediccion_E$ General_consensus$ mean<- prediccion_EC$ General_consensus$ mean<- prediccion_NE$ General_consensus$ mean
Comparación de los procedimientos
Cambios entre presente y futuro
Compararemos las áreas que permanecen adecuadas, se pierden o se ganan bajo cada procedimiento.
<- prediction_changes (current_predictions = mapa_actual,new_predictions = mapa_E,fitted_models = modelos_finales,predicted_to = "future" )<- prediction_changes (current_predictions = mapa_actual,new_predictions = mapa_EC,fitted_models = modelos_finales,predicted_to = "future" )<- prediction_changes (current_predictions = mapa_actual,new_predictions = mapa_NE,fitted_models = modelos_finales,predicted_to = "future" )
Cambios en la idoneidad ambiental
Guardar los resultados
dir.create ("Resultados" , showWarnings = FALSE )writeRaster (mapas, "Resultados/predicciones.tif" , overwrite = TRUE )writeRaster (cambios, "Resultados/cambios.tif" , overwrite = TRUE )
¿Qué identifica el MOP?
MOP compara las condiciones ambientales de calibración con las condiciones del escenario de transferencia .
Las condiciones de transferencia que quedan fuera de los rangos observados durante la calibración representan extrapolación estricta .
Calcular el MOP
<- single_mop (data = modelos_finales,new_variables = G_var,subset_variables = TRUE ,type = "detailed" ,na_in_range = FALSE ,calculate_distance = TRUE ,where_distance = "all" ,distance = "euclidean" ,scale = TRUE ,center = TRUE ,rescale_distance = TRUE ,progress_bar = FALSE
Condiciones no análogas
par (mfrow = c (1 , 3 ))plot (MOP_futuro$ mop_basic, main = "Condiciones no análogas" , col = c ("grey90" , "#D55E00" ), axes = FALSE )plot (MOP_futuro$ mop_simple, main = "Variables fuera del rango" , col = hcl.colors (8 , "YlOrRd" ), axes = FALSE )plot (MOP_futuro$ mop_distances, main = "Disimilitud ambiental" , col = hcl.colors (100 , "Inferno" ), axes = FALSE )
MOP básico: identifica celdas con al menos una variable fuera del rango de calibración.MOP simple: indica cuántas variables están fuera de sus rangos.Distancia ambiental: indica qué tan diferentes son las condiciones futuras.
Valores superiores al rango de calibración
plot ($ mop_detailed$ towards_high_end,main = paste (names (G_var), "por encima" ),col = c ("grey90" , "#D55E00" ),nc = 2 ,axes = FALSE )
Las celdas resaltadas presentan valores futuros superiores al máximo utilizado durante la calibración.
Valores inferiores al rango de calibración
plot ($ mop_detailed$ towards_low_end,main = paste (names (G_var), "por debajo" ),col = c ("grey90" , "#0072B2" ),nc = 2 ,axes = FALSE
Las celdas resaltadas presentan valores futuros inferiores al mínimo utilizado durante la calibración.
Interpretación del MOP
MOP identifica las regiones donde la transferencia ocurre bajo condiciones ambientales no representadas durante la calibración.
En esas regiones: - E: el modelo continúa sus respuestas fuera del rango; - EC: los valores se limitan al mínimo o máximo conocido; - NE: las celdas fuera del rango reciben idoneidad cero.
MOP indica dónde existe riesgo de extrapolación y la predicción debe interpretarse con mayor precaución.
¿Por qué algunas predicciones son mejores que otras?
¿Es posible mejorar las predicciones?
¿Por qué algunas predicciones son mejores que otras?
¿Es posible mejorar las predicciones?
