Objetivo

El objetivo de este proyecto es predecir el área de distribución de la Tenca (Mimus thenca), un ave endémica de Chile.

Para esto, crearé un Modelo de Distribución de Especies (SDM por sus siglas en inglés).

La idea es que en base al clima en que se ha encontrado la Tenca, se puede calcular la adecuabilidad de clima para esta especie en cada lugar, lo que funciona como un proxy para la probabilidad de encontrar esta especie a lo largo de todo Chile.

Existen muchos modelos posibles para proyectar la distribución de especies en base al clima, y el más utilizado es MaxEnt (Maximum Entropy) Phillips et. al, 2007. Para una guía práctica, ver Merow et. al, 2013

1. Librerías

Antes de empezar, cargo las librerías que necesitaré durante el resto del proyecto

library(raster) #para datos de clima
library(dismo) # para datos de ocurrencia y el modelo
library(rJava) # para la librería dismo
library(rgeos) # para la mascara de chile
library(grDevices) # para los colores del mapa

2. Input

Los datos necesarios para hacer un SDM son:

  1. Datos de clima del área a evaluar
  2. Datos de presencia de la especie

2.1. Datos de clima

Para determinar la distribución de la Tenca en Chile, usaré como predictores a las 5 variables climáticas que constringen la distribución geográfica de animales Araujo et. al, 2005. Sacaré esas variables de la pagina WorldClim, y luego restringiré los datos a la extensión de Chile.

presente <- getData("worldclim", var='bio', res=2.5) #obtiene todas las variables climáticas para todo el mundo
presente.a <- stack(presente, layers= c(1,10,11,12,18)) #RasterStack las 5 variables que me interesan
ext <-extent(-80,-65,-55.9,-17) #extensión chile
modelEnv <- crop(presente.a, ext) #datos de las 5 variables climáticas sólo para la extensión de chile

Además del clima, la topografía también modula la dispersión de especies, por lo que agregaré esa variable

topo <- raster("Topografia/gt30w100s10_dem/gt30w100s10.dem") # topografia del mundo
topo2 <- crop(topo, ext) # topografia de chile
topo3 <- resample(topo2, modelEnv) # en la misma resolución de pixeles que el clima
names(topo3) <- "Topo"  #nombre de la variable
modelEnv <- stack(modelEnv,topo3) # agrego la capa de topografía al RasterStack

2.2. Datos de presencia

En este caso, tenemos una tabla con las coordenadas en que se ha documentado presencia de la Tenca (columnas: Longitud y Latitud).

occs<-read.table("occ2019/mimus_thenca.txt",sep="\t",h=T)

2.3. Datos del Fondo (Background)

El modelo contrasta presencias contra localidades del “fondo” o “background”, por lo que ahora elijo puntos aleatorios del fondo. No se deben confundir con datos de ausencia, porque MaxEnt no ocupa ausencias. Para más detalles de los requerimientos y suposiciones de MaxEnt, ver Merow et. al, 2013

set.seed(0)
backGround <- randomPoints(modelEnv, 10000,p = occs,excludep = T) # 10000 puntos random del fondo, excluyendo puntos de presencia
colnames(backGround) <- c('Longitude', 'Latitude')

3. Datos Entrenamiento

3.1. Grupos de Entrenamiento y Prueba

Para poder evaluar cuán bueno es el modelo que crearé, usaré el 80% de los datos para entrenar el modelo y luego usaré el 20% restante para probar cuán bueno es el modelo (es decir, si efectivamente predice una alta adecuabilidad de hábitat en estos lugares donde sí hay presencia de la especie).

fold <- kfold(occs, k=5)  #dividir ocurrencias en 5 grupos 
modelTest <- occs[fold==2, ] # 20% de los datos serán para probar el modelo
modelTrain <- occs[fold!=2, ] # 80% de los datos serán para entrenar al modelo
trainPts <- rbind(modelTrain, backGround) # puntos de presencia y fondo para entrenar al modelo
presabs <- c(rep(1, nrow(modelTrain)), rep(0, nrow(backGround))) #vector con 1 = presencia y 0 = fondo

3.2. Clima de Entrenamiento

Ahora extraigo los valores de clima que tengo en modelEnv para cada punto con el que voy a entrenar mi modelo

trainVals <- as.data.frame(extract(modelEnv, trainPts))  #datos de clima de los puntos
para.m <- cbind(presabs, trainVals) #junto los valores de presencia/fondo para cada punto de entrenamiento
para.m <- na.omit(para.m) #omito los puntos donde alguna variable climática tiene NA
presabs.m <- para.m[,1] #recupero el vector de presencia/ausencia sin NA
trainVals.m <- para.m[,-1] # recupero el dataframe de datos de clima sin NA

4. MaxEnt

4.1. Crear Modelo

Con estos datos (puntos de presencia/fondo y valores climáticos de cada punto), ya puedo usar el método MaxEnt para hacer un modelo de distribución de la Tenca en Chile.

m <- maxent(p=presabs.m, x=trainVals.m) #objeto de class maxent
s <- predict(m,modelEnv) #predice la adecuabilidad de cada punto del mapa segun el modelo de maxent

4.2. Mapa Maxent

El objeto s indica la adecuabilidad del clima de cada pixel de Chile, por lo que lo podemos graficar:

load("CHL_adm0.RData") # crea el objeto gadm que es el poligono de chile
chile1<-crop(gadm,ext) #Máscara de Chile para que no aparezca un poco de argentina en el mapa
chi.s <- mask(s,chile1) #enmascara el raster con la capa de Chile, por lo que no se ve argentina
chile <- maps::map("world",regions = "Chile") #Poligono con la forma de chile para ver el contorno en el mapa
#Mapa
plot(chi.s,asp=1.4,cex.axis=0.8,legend=T,xlim=c(-90,-40),zlim=c(0,0.65),axis.args=list(cex.axis=0.8),
     legend.args=list(text="ROR",side=2),main="Distribución de Mimus thenca")
maps::map(chile,add=T, col="black")

ROR (Relative Occurrence Rate) describe la probabilidad relativa de que la celda esté contenida en la muestra de presencias. Esto significa que estos resultados se pueden interpretar como la adecuabilidad del hábitat (Habitat Suitability), no como la probabilidad de ocurrencia de una especie.

4.3 Evaluación del Modelo

Ahora podemos evaluar el modelo usando la funcion ‘evaluate’.

Los argumentos que necesita esta función son:
m = modelo maxent
p = presencias de testeo
a = fondo (Background) x = datos climáticos

e1 <- dismo::evaluate(m, p = modelTest, a=backGround, x=modelEnv)
e1
## class          : ModelEvaluation 
## n presences    : 2636 
## n absences     : 9982 
## AUC            : 0.9889864 
## cor            : 0.8593033 
## max TPR+TNR at : 0.3862824

El objeto e1 muestra el valor AUC (Área bajo la curva ROC). Este valor va entre 0 y 1. Un mayor valor indica un mejor modelo porque es más sensible (mayores verdaderos positivos) y específico (bajos falsos negativos). Éste modelo tiene un AUC = 0.99, por lo que es un modelo útil.

Bibliografía

Araújo, M.B. et, al. (2005) Validation of species–climate impact models under climate change. Global Change Biol. 11: 1504–1513.

Merow, C., Smith, M.J. & Silander, J.A. (2013) A practical guide to MaxEnt for modeling species’ distributions: what it does, and why inputs and settings matter. Ecography, 36: 1058–1069.

Phillips, S. et. al. (2006) Maximum entropy modeling of species geographic distributions. Ecol. Model. 190: 231–259.