En las dos últimas décadas, las variables bioclimáticas han sido utilizadas como una herramienta importante para entender la distribución de especies y priorizar áreas para la conservación de las especies a través del Modelado de Nicho Ecológico (MNE).
Por lo cual las variables bioclimáticas han cobrado relevancia como variables predictivas de la distribución de las especies.
Sin embargo, los conjuntos de datos interpolados de las variables bioclimáticas causan un sobreajuste de los modelos, debido principalmente a la multicolinealidad o redundancia dentro de las variables.
• Realizar la selección de variables empleando el VIF, como una herramienta para su uso en modelos de distribución de especies.
• Llevar a cabo un ánalisis de componentes principales de variables descriptoras de la distribución de especies.
• Mostrar la bondad de R, en el análisis de grandes volumenes de información; integrando procesos de obtención y manipulación de datos, análisis estadístico y SIG´s.
• Emplear el modelo de “nicho climático” o “nicho grinelliano” para tener un proxy de la distribución de cada especie de la subsección Oocarpae partir del espacio ambiental
• Crear un fichero html usando markdown para su publicación on-line, como estratégia para compartir y retroalimentar código de R aplicado a un área en particular.
• Poner en práctica conocimientos adquiridos, así como el uso de herramientas de cómputo, abordados en el curso.
# Remoción de objetos del ambiente de trabajo (Limpieza)
rm(list=ls())
# Directorio de trabajo
setwd("I:/Nueva carpeta")# Carga de paquetes a emplear
#
library(sp) # Manipulación de datos espaciales
library(rgdal) # Operaciones con datos Geoespaciales
library(raster) # Creación de archivos raster
library(dismo) # Modelos de distribución de especies
library(maptools) # Combinar datos espaciales
library(usdm) #
library(ade4)
library(ape) # Análisis de filogenética y evolución
library(ecospat)
library(gbm)
library(ecospat) # Métodos en ecología espacialVariables bioclimaticas de Worldclim ajustadas para México por Cuervo et al, 2013. Las cuales se obtuvieron de la interpolación de datos meteorológicos del periodo de 1910 - 2009. Datos en formatos matriciales, tiff y rst para desplegar directamente en SIG. La resolución espacial es de 30 arcos de segundo (≈1 km),con sistema de referencia latitud/longitud (no proyectadas) y datum WGS84.
# Remoción de objetos del ambiente de trabajo (Limpieza)
# rm(list=ls())# Descarga remota de dataset variables bioclimáticas ajustadas para México
# download.file("http://idrisi.uaemex.mx/images/Clima/BioClim(Tiff).rar", "Nueva carpeta", method = "auto")
# Extracción de archivos de variables bioclimáticas
# unzip(BioClim(Tiff))Cuervo-Robayo A. P., O. Téllez-Valdés, M. Gómez, C. Venegas-Barrera, J. Manjarrez y E. Martínez-Meyer. (2013). An update of high-resolution monthly climate surfaces for Mexico. International Journal of Climatology. Doi: 10.1002/joc.3848. email: ancuervo@gmail.com link: http://idrisi.uaemex.mx/distribucion/superficies-climaticas-para-mexico
# BIO1 Temperatura media anual
# BIO2 Rango de temperaturas diurnas
# BIO3 Isotermalidad (BIO2/BIO7) (* 100)
# BIO4 Estacionalidad en la temperatura (desviación estándar * 100)
# BIO5 Temperatura máxima del mes más cálido
# BIO6 Temperatura mínima del mes más frío
# BIO7 Rango anual de temperatura (BIO5-BIO6)
# BIO8 Temperatura media del trimestre más lluvioso
# BIO9 Temperatura media del trimestre más seco
# BIO10 Temperatura media del trimestre más cálido
# BIO11 Temperatura media del trimestre más frío
# BIO12 Precipitación anual
# BIO13 Precipitación del mes más lluvioso
# BIO14 Precipitación del mes más seco
# BIO15 Estacionalidad en la precipitación (coeficiente de variación)
# BIO16 Precipitación del trimestre más lluvioso
# BIO17 Precipitación del trimestre más seco
# BIO18 Precipitación del trimestre más cálido
# BIO19 Precipitación del trimestre más fríoConvierte un archivo ASCII que representa datos ráster en un dataset ráster
bio1 <- raster("J:/Research_2019/Bioclim_2019/ASCII_BioclimaticVariables/bio1.asc")
bio2 <- raster("J:/Research_2019/Bioclim_2019/ASCII_BioclimaticVariables/bio2.asc")
bio3 <- raster("J:/Research_2019/Bioclim_2019/ASCII_BioclimaticVariables/bio3.asc")
bio4 <- raster("J:/Research_2019/Bioclim_2019/ASCII_BioclimaticVariables/bio4.asc")
bio5 <- raster("J:/Research_2019/Bioclim_2019/ASCII_BioclimaticVariables/bio5.asc")
bio6 <- raster("J:/Research_2019/Bioclim_2019/ASCII_BioclimaticVariables/bio6.asc")
bio7 <- raster("J:/Research_2019/Bioclim_2019/ASCII_BioclimaticVariables/bio7.asc")
bio8 <- raster("J:/Research_2019/Bioclim_2019/ASCII_BioclimaticVariables/bio8.asc")
bio9 <- raster("J:/Research_2019/Bioclim_2019/ASCII_BioclimaticVariables/bio9.asc")
bio10 <- raster("J:/Research_2019/Bioclim_2019/ASCII_BioclimaticVariables/bio10.asc")
bio11 <- raster("J:/Research_2019/Bioclim_2019/ASCII_BioclimaticVariables/bio11.asc")
bio12 <- raster("J:/Research_2019/Bioclim_2019/ASCII_BioclimaticVariables/bio12.asc")
bio13 <- raster("J:/Research_2019/Bioclim_2019/ASCII_BioclimaticVariables/bio13.asc")
bio14 <- raster("J:/Research_2019/Bioclim_2019/ASCII_BioclimaticVariables/bio14.asc")
bio15 <- raster("J:/Research_2019/Bioclim_2019/ASCII_BioclimaticVariables/bio15.asc")
bio16 <- raster("J:/Research_2019/Bioclim_2019/ASCII_BioclimaticVariables/bio16.asc")
bio17 <- raster("J:/Research_2019/Bioclim_2019/ASCII_BioclimaticVariables/bio17.asc")
bio18 <- raster("J:/Research_2019/Bioclim_2019/ASCII_BioclimaticVariables/bio18.asc")
bio19 <- raster("J:/Research_2019/Bioclim_2019/ASCII_BioclimaticVariables/bio19.asc")# Despliegue de mapas raster
# > par(mfrow=c(4,5))
# > plot(bio1, main="Temperatura media anual")
# > plot(bio2, main="Rango de temperatura media diurna
# (Tmax-Tmin)")
# > plot(bio3, main="Isotermalidad
# (bio1/bio7 * 100)")
# > plot(bio4, main="Estacionalidad de T
# (DS * 100)")
# > plot(bio5, main="T máxima del mes más caliente")
# > plot(bio6, main="T mínima del mes más frío")
# > plot(bio7, main="Rango de T anual (bio5-bio6)")
# > plot(bio8, main="T media del trimestre más húmedo")
# > plot(bio9, main="T media del trimestre más seco")
# > plot(bio10, main="T media del trimestre más frío")
# > plot(bio11, main="T media del trimestre más caliente")
# > plot(bio12, main="Precipitación total anual")
# > plot(bio13, main="Precipitación del mes más húmedo")
# > plot(bio14, main="Precipitación del mes más seco")
# > plot(bio15, main="Estacionalidad de la precipitación")
# > plot(bio16, main="Precipitación trimestre más húmedo")
# > plot(bio17, main="Precipitación trimestre más seco")
# > plot(bio18, main="Precipitación trimestre más caliente")
# > plot(bio19, main="Precipitación trimestre más frío")Apilado de raster para extracción de la información bioclimatica de las localidades
# > mexico.stk <- stack(bio1, bio2, bio3, bio4, bio5, bio6,
# > bio7, bio8, bio9, bio10, bio11, bio12,
# > bio13, bio14, bio15, bio16, bio17, bio18, bio19)# Eliminar raster individuales
# rm(bio1,...bio19)Aproximadamente 40 % de las especies de pino se encuentran en México. De las cuales el grupo de pinos oocarpa presenta el mayor rango de distribución norte-sur de todas las especies de pinos presentes en México y Centro América, extendiéndose desde el noroeste de México (e.g. sur de Sonora), hacia el sur a través de México y Guatemala, hasta Belice, El Salvador, Honduras y el noroeste de Nicaragua (Richardson, 1998; Perry, 1991).
Se obtuvieron registros de presencia de la especies Pinus oocarpa de las colectas realizadas por CAMCORE.
Dvorak, W.S. et al (2000) The evolutionary history of Mesoamerican Oocarpae. In: Conservation & Testing of Tropical & subtropical Forest Tree Species by CAMCORE Cooperative. 1-11.
# Pinus oocarpa
sp_1 <- read.csv("I:/Nueva carpeta/Localidades_P_oocarpa.csv",
na.strings = 0)
# View(sp_1)
sp1_points <- sp_1[,c("longitud","latitud")]
# plot(sp1_points[,1],sp1_points[,2],col="blue",pch=19)# Explorar el tamaño del archivo raster, resolución, clase, extensión y proyección
bbox(bio1); ncol(bio1); nrow(bio1); res(bio1); class(bio1); extent(bio1); projection(bio1)## min max
## s1 -122 -79
## s2 13 34## [1] 5160## [1] 2520## [1] 0.008333333 0.008333333## [1] "RasterLayer"
## attr(,"package")
## [1] "raster"## class : Extent
## xmin : -122
## xmax : -79
## ymin : 13
## ymax : 34## [1] "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"Extensión del raster: Latitud 13 a 34 y Longitud de -79 a -122 grados. Resolución = 30 arc segundo = 30.9 x 30 m = 927 m aprox. 1 km x 1 km
# Definir extensión del área de trabajo
map.ext <- extent(-122, -79, 13, 34)# Extraer y graficar capas del objeto stack
# plot(mexico.stk$bio1, col=rainbow(100, start=.0, end=.8), ext=map.ext)# sp1_points_mat <- as.matrix(sp1_points)
# sp1_dups <- duplicated(sp2_points)
# sp1_dups_row <- which(sp2_dups==TRUE)
# sp1_points_mat <- sp2_points[-sp2_dups_row,]
# sp1_points_mat <- na.omit(sp2_points_mat)library(ecospat)
# > ecospat.cor.plot(mexico.stk[ ,1:19])
# > vifstep(mexico.stk[[1:19]])
# > vifcor(mexico.stk, th=.7)La mayoría de las técnicas estadísticas tendrán dificultades para ajustarse con éxito a un modelo estable, si las variables predictoras están altamente correlacionadas entre sí. (Guisan et al, 2017)
Existen métodos para evitar lo anterior, una manera es considerar sólo algunas variables predictoras, lo que hace resulta en modelos con “parsimonia” y menor sobre/sub estimación.
¿Cómo realizar la remoción de variables que estén fuertemente correlacionadas?
En diversos artículos de modelos de distribución han tomado como umbral valores de r=0.8 y r=0.7 (Menard, 2002; Green, 1979)
Métodos para minimizar la redundancia de las variables bioclimáticas predictoras incluyen:
El VIF estima la severidad del efecto de la multicolinealidad, midiendo el grado en que la varianza en una regresión aumenta debido a la colinealidad, en comparación con el uso de variables no correlacionadas.
Función VIF() vegan package y usdm package en R
En este caso valores de VIF cercanos a 10 son considerados como críticos en la correlación de variables múltiples. (Guisan et al, 2017)
# Factor de inflación de varianza
# >library(usdm)
# Stepwise eliminacion de variables con inflación, con un valor default de VIF=10
# > vifstep(env.variables)
# Identificar variables con presencia de colinearidad r=0.7
# > vifcor(env.variables, th=.7)Una de las aplicaciones de PCA es la reducción de dimensionalidad (variables), perdiendo la menor cantidad de información (varianza) posible.
Útil cuando contamos con un gran número de variables cuantitativas posiblemente correlacionadas (indicativo de existencia de información redundante).
PCA permite reducirlas a un número menor de variables transformadas (componentes principales) que expliquen gran parte de la variabilidad en los datos.
# RASTER PCA
library(RStoolbox)
# > pca<-rasterPCA(mexico.stk , nSamples = NULL, nComp = nlayers(mexico.stk), spca = FALSE, maskCheck = TRUE)
# > summary(pca$model)
# > knitr::kable(round(pca$model$loadings[, 1:4],3))
# > plot(pca$map,1, legend=F)
# > plot(pca$map,2, legend=F)
# > plot(pca$map,3, legend=F)
# > plot(pca$map,4, legend=F)
# > pca$model$loadings
# Convertir a raster
# > pca1<-raster(pca$map,1)
# > pca2<-raster(pca$map,2)
# > pca3<-raster(pca$map,3)
# > pca4<-raster(pca$map,4)
# > pcas<-stack(pca1,pca2, pca3, pca4)
# Imprimir raster de componentes principales
# > ggRGB(pca$map,1,2,3,4, stretch="lin", q=0)
# > if(require(gridExtra)){
# > plots <- lapply(1:4, function(x) ggR(pca$map, x, geom_raster = TRUE))
# > grid.arrange(plots[[1]],plots[[2]], plots[[3]], plots[[4]], ncol=2)}# > mex.env <- crop(mexico.stk, mex.extent)
# > mexico.stk3<- stack (bio2, bio3, bio8, bio9, bio10 , bio15, bio16, bio18, bio19)
# > mex.env3 <- crop(mexico.stk3, mex.extent)
# > mex.env4 <- crop(pcas, mex.extent)# > sp1_bioclim <- bioclim(mex.env, sp1_points_mat)
# > plot(sp1_bioclim, a=2, b=6, p=0.85)# sp1.map<- predict (mex.env, sp1_bioclim)group <- kfold(sp1_points, 5)
pres_train <- sp1_points[group != 1, ]
pres_test <- sp1_points[group == 1, ]
#backg <- randomPoints(mex.env, n=1000, ext=mex.extent,
#extf = 1.25)
#colnames(backg) = c('longitud', 'latitud')
#group4 <- kfold(backg, 5)
#backg_train <- backg4[group != 1, ]
#backg_test <- backg[group == 1, ]
#eval_test <- evaluate(pres_test, backg_test, sp1_bioclim, mex.env)
#eval_testMapas de distribución: Obtención de mapas de distribución geográfica para cada especie de la subsección Oocarpae, a partir de registros de presencia del dataset de CAMCORE-NCSU y de variables bioclimáticas para México.
# > plot (sp1.map, main='Pinus oocarpa: áreas favorables
# (19 variables) ModelEvaluation AUC:0.9625')
# Agregar los puntos "reales" de la especie
# > points(sp1_points_mat[,1],sp1_points_mat[,2],cex=0.4)
# > plot (sp1.map2, main='Pinus oocarpa: áreas favorables
# (VIF & Corr variables)')
# Agregar los puntos "reales" de la especie
# > points(sp1_points_mat[,1],sp1_points_mat[,2],cex=0.4)
# > plot (sp1.map3, main='Pinus oocarpa: áreas favorables
# (VIF variables) ModelEvaluation AUC:0.9492')
# Agregar los puntos "reales" de la especie
# > points(sp1_points_mat[,1],sp1_points_mat[,2],cex=0.4)
# > plot (sp1.map4, main='Pinus oocarpa: áreas favorables
# (PCAs) ModelEvaluation AUC:0.9085')
# Agregar los puntos "reales" de la especie
# > points(sp1_points_mat[,1],sp1_points_mat[,2],cex=0.4)Establecer un “threshold” para cortar la predicción y generar un mapa binario
# > threshold <- e@t[which.max(e@TPR + e@TNR)]
# > plot(sp1.map > threshold, main='Pinus oocarpa: presencia/ausencia')
# Agregar los puntos "reales" de la especie
# > points(sp1_points_mat[,1],sp1_points_mat[,2],cex=0.4)El uso de variables seleccionadas en función del VIF, si bien presentó el segundo mejor desempeño (AUC=0.94), se ajusto mejor en el espacio geográfico a la distribución real de la especie.
Trabajar sólo con variables no redundantes no siempre da buenos resultados, ya que algunas de éstas variables pueden actuar como un buen descriptor ecológico.
El modelo final puede ser seleccionado en base al más alto valor de AUC e incorporando un menor número de variables correlacionadas (r>0.8, r2>0.8).
Resulta útil seguir el consejo del experto de la especie/género/familia en la elección de variables, lo que puede contribuir a incorporar variables limitantes y tienen un significado biológico en la distribución de la especie.