1 Comparación de técnicas para la obtención de modelos de distribución de especies con RVII JORNADAS DE USUARIOS DE R Salamanca 2015 Jennifer Morales Barbero* Dolores Ferrer Castán VII JORNADAS DE USUARIOS DE R Salamanca 2015 Jennifer Morales Barbero* Dolores Ferrer Castán

2 Introducción Hutchinson (1957) definía el nicho ecológico de una especie como un hipervolumen de n-dimensiones en el que se dan las condiciones ambientales necesarias para que la especie puede sobrevivir. Áreas con una combinación similar de valores de esas mismas variables permitirían la presencia de la especie, generando su distribución potencial. Los modelos de distribución de especies son una representación parcial de la realidad que trata de generar mapas que representen la idoneidad de un espacio para una especie. Nicho realizado Nicho potencial

3 Esquema Modelado de distribución de especiesIntroducción Esquema Modelado de distribución de especies Datos de entrada Modelado de nicho Predicción de la distribución potencial de la especie Variable 3 Variable 2 Variable 1 Evaluación del modelo

4 Introducción Asunciones:Las variables ambientales determinan el nicho ecológico. La especie se encuentra en equilibrio o pseudoequilibrio con el ambiente. Incertidumbres: Calidad de capas y datos disponibles. Exactitud, complejidad y naturaleza del modelo. No consideran interacciones bióticas, barreras geográficas, capacidad de dispersión, etc. Aplicaciones: Cartografía automática de especies. Análisis de distribución de especies. Localizaciones de poblaciones desconocidas. Predicción de invasiones biológicas. Reconstrucción de distribuciones del pasado. Predicción de efectos del cambio climático. …

5 Variables ambientalesDatos de entrada Variables ambientales library(raster) library(maps) library(mapdata) Importar los datos de las variables predictoras. WorldClim – Global Climate Data Confirmar extensión, resolución, formato y proyección. bio1<-raster(“C:/.../bio1.asc") # extensión Península Ibérica, formato .asc ... predictores <- stack(bio1, bio3, bio5, bio7, bio12, bio14, bio15) class : RasterStack dimensions : 240, 336, 80640, 7 (nrow, ncol, ncell, nlayers) resolution : , (x, y) extent : -10, 4, 35, 45 (xmin, xmax, ymin, ymax) coord. ref. : +proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs names : bio1, bio5, bio7, bio14, bio15, bio3, bio12 min values : -7, 113, 123, 0, 11, 26, 219 max values : 193, 364, 335, 105, 86, 50, 1757 predictores_6k <- stack(bio1, bio3, bio5, bio7, bio12, bio14, bio15) predictores_2070<- stack(bio1, bio3, bio5, bio7, bio12, bio14, bio15) # extensión Península Ibérica, formato.asc, renombrado a bio1,...

6 Variables ambientalesDatos de entrada Variables ambientales library(raster) library(maps) library(mapdata) Environmental data WorldClim – Global Climate Data bio1<-raster(“C:/.../bio1.asc") # extension Península IbéricaI, format .asc ... predictors <- stack(bio1, bio3, bio5, bio7, bio12, bio14, bio15) class : RasterStack dimensions : 240, 336, 80640, 7 (nrow, ncol, ncell, nlayers) resolution : , (x, y) extent : -10, 4, 35, 45 (xmin, xmax, ymin, ymax) coord. ref. : +proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs names : bio1, bio5, bio7, bio14, bio15, bio3, bio12 min values : -7, , 123, 0, , , 219 max values : 193, 364, 335, , , , 1757 predictors_6k <- stack(bio1, bio3, bio5, bio7, bio12, bio14, bio15) # extensión PI, format.asc, renombrado a bio1,...

7 Variable respuesta (Presencias)Datos de entrada Variable respuesta (Presencias) Importar datos de ocurrencia de la especie Global Biodiversity Information Facility (http://www.gbif.org) Limpieza de datos Eliminar registros duplicados sp.occ <-read.csv("C:/…/Sorex_granarius.csv", header=T, sep=",") # formato.csv. 2 columnas (Longitude/Latitude). 2 decimales. coordinates(sp.occ)<-c("longitude","latitude") presclim <- extract(predictores, sp.occ, method='bilinear', buffer=NULL, fun=NULL, df=TRUE) #Extraer valores del raster. pres=rep(1, nrow(presclim)) presclim = data.frame(presclim, pres) presencias<-cbind(presclim,sp.occ)

8 Variable respuesta (“background points”)Datos de entrada Variable respuesta (“background points”) library(dismo) #cf Dolores Ferrer Castán abs <-read.csv("C:/…/backgroundpoints.csv", header=T) # 500 puntos background, formato .csv; latitud y longitud, redondeado a 2 dígitos; sin duplicados coordinates(abs)<- c("longitude","latitude") absclim <- extract(predictores, abs, method='bilinear', buffer=NULL, fun=NULL, df=TRUE) # extraer valores del raster pres = rep(0, nrow(absclim)) absclim = data.frame(absclim, pres) ausencias <- cbind(absclim, background) presabs <- rbind(presencias, ausencias) presabs = data.frame(presabs) presabs[!duplicated(presabs[,8:9]),] # [,columna longitude:latitude] ID bio1 bio2 bio7 bio12 bio15 pres longitude latitude Estudio de correlación entre variables cor.test(variable1, variable2, alternative="two.sided", method="pearson") https://ecologicaconciencia.wordpress.com/

9 Variable respuesta (“background points”)Datos de entrada Variable respuesta (“background points”)

10 Modelado de nicho Tipos de algoritmoMétodos basados en árboles de decisión: Random Forest (RF; Breiman, 2001). Selecciona puntos al azar (muestreo con reemplazo) para crear diferentes sets de datos Al crear los árboles, las variables se eligen en cada nodo del árbol Crea un árbol de decisión con cada set de datos, obteniendo diferentes árboles (cada set contiene diferentes puntos y diferentes variables) Predice los nuevos datos usando el "voto mayoritario", donde clasificará como "positivo" si la mayoría de los árboles predicen la observación como positiva Datos de presencia y “ausencia” 500 “background points” 500 árboles. Modelado de nicho Métodos de regresión: Modelos Aditivos Generalizados (GAMs; Hastie y Tibshirani, 1990; Yee y Mitchell, 1991) Modelos Lineales Generalizados (GLM; McCullagh y Nelder, 1989)) Distribución binomial de los errores y vinculación logística. Selección de variables paso a paso hacia delante Datos de presencia y “ausencia”. 500 “background points” Tipos de algoritmo

11 Modelado de nicho GAMs Probabilidad de ocurrencia de la especieRelación entre la presencia de Sorex granarius con las variables ambientales más relevantes. Las líneas de regresión han sido generadas mediante modelos aditivos generalizados (GAMs) del paquete “gam”. BIO1 BIO5 BIO7 BIO12 BIO14 BIO15 Probabilidad de ocurrencia de la especie P < 0,0001 P < 0,0001 P < 0,0001 P < 0,01 P < 0,0001 P < 0,0001

12 Modelado de nicho Evaluación Matriz de confusión Curva ROC(Receiver Operating Characteristic) AUC (Area under the curve) (0.5-1) [0.5, 0.6): Test malo. [0.6, 0.75): Test regular. [0.75, 0.9): Test bueno. [0.9, 0.97): Test muy bueno. [0.97, 1): Test excelente. Evaluación Presencia real Ausencia real Presencia predicha A Verdadero positivo B Falso negativo Error Comisión (sobrepredicción) Ausencia predicha C Error Omisión (subpredicción) D Verdadero negativo A/(A+C) SENSIBILIDAD (Fracción de verdaderos positivos) C/(A+C) TASA OMISIÓN (Fracción de falsos positivos) D/(D+B) ESPECIFICIDAD (Fracción de verdaderos negativos)

13 Métodos de modelado y EvaluaciónModelado de nicho GLMs (paso a paso hacia adelante) modelo <- glm (pres ~ I(bio14^2)+bio1+I(bio7^2)+I(bio7^3)+bio5+bio15+I(bio15^2)+I(bio15^3), binomial, data=presabs) presencias <-read.table("C:/…/sorex_granarius.txt", header=T, dec=",")# 2 columnas (longitude;latitude) ausencias <-read.table("C:/…/backgroundpoints.txt", header=T, dec=",") coordinates(presencias)<-c(“longitude", “latitude") coordinates(ausencias)<-c(“longitude", “latitude") e <- evaluate(presencias, ausencias, modelo, predictores) e class : ModelEvaluation n presences : 172 n absences : 469 AUC : cor : max TPR+TNR at : plot(e, 'ROC', col='blue', cex=0.1) Métodos de modelado y Evaluación

14 Métodos de modelado y EvaluaciónModelado de nicho Random Forest library(randomForest) modelo <- pres ~ bio1 + bio3 + bio5 + bio7 + bio12 + bio14 + bio15 #variables no correlacionadas RF <- randomForest(modelo, data=presabs) e <- evaluate(presencias, ausencias, RF, predictores) e class : ModelEvaluation n presences : 172 n absences : 469 AUC : cor : max TPR+TNR at : plot(e, 'ROC', col=‘red', cex=0.1) Métodos de modelado y Evaluación

15 Modelado de nicho Predicción de Modelosnumerador <- exp( e+01-( e-03*bio14*bio14)-( e-01*bio1)-( e-03*bio7*bio7)+( e-06*bio7*bio7*bio7)+( e-01*bio5)-( e+00*bio15)+( e-02*bio15*bio15)-( e-04*bio15*bio15*bio15)) denominador <- 1+numerador cociente <- numerador/denominador plot(cociente) Modelado de nicho GLMs (paso a paso hacia adelante) modelo <- glm (pres ~ I(bio14^2)+bio1+I(bio7^2)+I(bio7^3)+bio5+bio15+I(bio15^2)+I(bio15^3), binomial, data=presabs) predictores <- stack(bio1, bio3, bio5, bio7, bio12, bio14, bio15) summary(modelo) Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 3.433e e ** I(bio14^2) e e e-12 *** bio e e e-11 *** I(bio7^2) e e e-11 *** I(bio7^3) 2.042e e e-08 *** bio e e e-07 *** bio e e *** I(bio15^2) 7.347e e *** I(bio15^3) e e e-05 *** --- Null deviance: on 640 degrees of freedom Residual deviance: on 632 degrees of freedom AIC: ( )*100/745.6 [1] # %varianza explicado Predicción de Modelos

16 Modelado de nicho Predicción de Modelos Random Forestp<-predict(predictores, RF) plot(p) Modelado de nicho Random Forest modelo <- pres ~ bio1 + bio3 + bio5 + bio7 + bio12 + bio14 + bio15 predictores <- stack(bio1, bio3, bio5, bio7, bio12, bio14, bio15) RF <- randomForest(modelo, data=presabs) Coefficients: (Intercept) bio bio5 bio7 bio14 bio15 Degrees of Freedom: 640 Total (i.e. Null); 635 Residual Null Deviance: Residual Deviance: AIC: RF Call: randomForest(formula = model, data = presabs) Type of random forest: regression Number of trees: 500 No. of variables tried at each split: 2 Mean of squared residuals: % Var explained: 65.47 Predicción de Modelos

17 Proyección de modelos en espacio y tiempoModelado de nicho GLMs (modelo hacia adelante) modelo <- glm (pres ~ I(bio14^2)+bio1+I(bio7^2)+I(bio7^3)+bio5+bio15+I(bio15^2)+I(bio15^3), binomial, data=presabs) predictores_6k <- stack(bio1, bio3, bio5, bio7, bio12, bio14, bio15) summary(modelo) Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 3.433e e ** I(bio14^2) e e e-12 *** Bio e e e-11 *** I(bio7^2) e e e-11 *** I(bio7^3) e e e-08 *** bio e e e-07 *** bio e e *** I(bio15^2) e e *** I(bio15^3) e e e-05 *** Proyección de modelos en espacio y tiempo #Opción 1: modelo con coeficientes de presente. numerador <- exp( e+01-( e-03*bio14*bio14)-( e-01*bio1)-( e-03*bio7*bio7)+( e-06*bio7*bio7*bio7)+( e-01*bio5)-( e+00*bio15)+( e-02*bio15*bio15)-( e-04*bio15*bio15*bio15)) denominador <- 1+numerador cociente <- numerador/denominador plot(cociente) #Opción 2: función predict pred <- predict(predictores_6k, modelo, type="response") plot(pred)

18 Proyección de modelos en espacio y tiempoModelado de nicho Random Forest modelo <- pres ~ bio1+bio3+bio5+bio7+bio12+bio14+bio15 RF <- randomForest(modelo, data=presabs) predictors_6k <-stack(bio1, bio3, bio5, bio7, bio12, bio14, bio15) pred <- predict(predictores_6k, RF, type="response") plot(pred) Proyección de modelos en espacio y tiempo Modelos consenso: models<-(stack(pred1, pred2, predN)) plot(mean(models))

19 (paso a paso hacia delante)Modelado de nicho Binario e=evaluate(presencias, ausencias, modelo, predictores) pr <-predict(predictores, modelo) tr<-threshold(e, 'spec_sens') # umbral en el cual la suma de sensibibidad y especificidad es mayor. bin<-pr > tr plot(bin) GLMs (paso a paso hacia delante) Random Forest

20 (paso a paso hacia adelante)Modelado de nicho Resumen Resultados Medio Holoceno Presente (paso a paso hacia adelante) GLMs AUC=0.943 56.01 %varianza explicado AUC=0.994 65.47 %varianza explicado Random Forest

21 ¡¡MUCHAS GRACIAS por su atención!!