Introducción

En este módulo se busca explorar las herramientas de R para el manejo de datos biológicos.Se va a construir un modelo de distribución de especie (Species Distribution Model [SDM]) para el oso de anteojos Tremarctos ornatus. Se van a responder dos preguntas:

¿Cómo es el nicho ecológico de T. ornatus?
¿Cómo está distribuida la diversidad (riqueza de especies) de los mamñiferos terrestres carnívoros en Colombia y en el mundo?

Modelo de distribución de especies

Formulación que relaciona variables bióticas y abióicas de un paisaje con los sitios de tal paisaje donde ha sido avistada una especie. A través de este modelo se busca:

Explorar las variables que son explicativas de de la presencia de la especie.
Predecir dónde puede estar la especie.

Entonces lo que se necesita para hacer el SDM del oso de anteojos es la información de su presencia, regristros georreferenciados (latitid, longitud y altitud), e información del ambeinte donde se da el avistamiendo, capas ambientales. Con esto se busca medir a probabilidad de ocurrencia del oso en el paisaje.

Internacionalmente existen varias plataformas que organizan y almacenan registros georreferenciados, la más popular y la que usaremos en este módulo es GBIF (Global Bioinformation Informataion Facility); para Colombia existe el SIB, Sistema de Inforamción sobre Biodiversidad. Pero por convenio, todo lo que esté en el SIB está en el GBIF.

Datos georreferenciados

Para obtener datos de GBIF se puede ir directamente a su página o se puede usar el paquete de R rgbif (diseñado por hattps://ropensci.org/) o con el paquete dismo mediante la función gbif(). El primer paquete se comunica con GBIF. por ejemplo, para traer la infroamción del oso de anteojos se puede:

library(rgbif)
oso = occ_search(
  scientificName = "Tremarctos ornatus", 
  limit = 2000, 
  country = "CO"
  hasCoodinate = TRUE,
  hasGeospatialissue = FALSE,
  fields = "minimal",
  )$data
writecsv2(oso, "Oso,csv")

Otra opción es muestrar uno mismo de manera que se tenga la posibilidad de determinar las ocurrencias y las ausencias.

Datos de información abiental terrestre

Esta infomación suele venir en formato ráster (mediante pixeles). El depósito internación más completo a nivel global es el de BIOCLIM. En esta plataforma se puede encontrar variables como isotermalidad, temperatura media anual, precipitación anual, en totoal hay 19 capas. También se puede encontrar información para el futuro según modelos climáticos. Las 19 varaibles son:

BIO1 = Annual Mean Temperature
BIO2 = Mean Diurnal Range (Mean of monthly (max temp - min temp))
BIO3 = Isothermality (BIO2/BIO7) (×100)
BIO4 = Temperature Seasonality (standard deviation ×100)
BIO5 = Max Temperature of Warmest Month
BIO6 = Min Temperature of Coldest Month
BIO7 = Temperature Annual Range (BIO5-BIO6)
BIO8 = Mean Temperature of Wettest Quarter
BIO9 = Mean Temperature of Driest Quarter
BIO10 = Mean Temperature of Warmest Quarter
BIO11 = Mean Temperature of Coldest Quarter
BIO12 = Annual Precipitation
BIO13 = Precipitation of Wettest Month
BIO14 = Precipitation of Driest Month
BIO15 = Precipitation Seasonality (Coefficient of Variation)
BIO16 = Precipitation of Wettest Quarter
BIO17 = Precipitation of Driest Quarter
BIO18 = Precipitation of Warmest Quarter
BIO19 = Precipitation of Coldest Quarter

Unión de datos georreferenciados y ambientales

Para unir los datos en un modelo de dsitribución de especies usamos el paquete sdm; empero si queremos utilizar una GUI podemos ir a https://wallaceecomod.github.io/. Una vez instalado sdm, se carga la librería y se instala todos los paquetes asociados con sdm::installA().

library(sdm)
installAll()

Evalución del modelo

Una vez se tienen el modelo lo podemos evaluar particionando los datos con datos de entrenamiento y datos de prueba. Se crea un modelo con los datos de entrenamiento y con este se intenta predecir los datos de prueba, segun qué tan correctamente sean las prediciones se evalua el desempeño del modelo. Otra manera es usar validación cruzada (corss validation), recurrentemente se apartan unos datos y se predicen con los otros datos.

Normalmente lo que surge es una matriz de confusión:

True positives (TP)	False positives (FP)
True negatives TN	False negatives (FN)

A partir de la cual se obtiene la tasa de verdaderos positivos, \(TPR\), y la tasa de falsos positivos, \(FPR\).

\[ TPR = \frac{TP}{TP+FN} \]

\[ FPR = \frac{FP}{FP+TN} \]

Positivo representa presencias y negativos representan ausencias. Con dichas tasas se puede crear un umbral de presencias y ausencias para la especie, debajo o arriba de cada umbral hay una probabilida y una tasa de que la especie esté o no. La forma de visualizar todos los posibles umbrales de un modelo es mediante el receiver operating characteristic curve (ROC), en que se grafican las tasas de falsos positivos por las tasas de verdaderos positivos. Con esta curva se puede saber cuál sería el mejor umbral, el que maximice los verdaderos positivos.

Glosario

Sistema de referencia de coordenadas

Se usa para visualizar el geoide terrestre de manera plana. Esto implica una pryección a dos dimensiones, lo cuel es una distorsión que hay que tener en cuenta. Detrás de la proyección hay un modelo matemático que se resumen en el datum. El más común es el WGS84.También se usa UTM. En particular, en Colombia se usa el Magna-Sirgas.

Vif (Variable Inflaction Factor)

Procedimiento para medir la colinearidad (cuando variables independientes y predictorias [ambiente] están correlacionadas) en un conjunto de datos de variables predictorias. Se escogen las varaibles predictorias (ambientales) y se toma una para hcer que todas las otras variables la predigan, si la predicen bien puede tratarse de una variable redundante.

Seudoausencias

La variable que se modela es binomial (presencia-ausencia). Varias algoritmos, en particular los de regresión, requieren ausencia, por lo que se genera un número aleatorio de ausencias.

Algoritmos

Analizan modelos de distribución de especies teneindo en cuenta el esquema BAM, que define dónde peude estar una especie: en la intersección entre moviento, varaibles bióticas y variables abióticas.

De envoltura (Envelopes)

Crean una especie de objeto multidimensional donde cada dimensión es una variable ambiental con un mínimo y un máximo. Se prioduce una región entre los mínimos y máximos donde se encuentra la especie.

bioclim
bioclim.dismo
domain.dismo
mahal.dismo

De regresión

Modelos lineares generalizados
Modelos geenralizados aditivos
Análisis flexible discriminante etc.

De aprendizaje autómata (machine learning)

Existen procedimeintos iterativos que en función de los datos buscan predecir las presencias y ausencias.

Random forest
Boosted regression tree
support-vector machines etc.

Práctica: Tremarctos ornatus en el mundo

Comenzamos cargando los paquetes. El paquete sdm llama el paquete sp y el paquete dismo carga a su vez el paquete raster.

library(sdm)

## Loading required package: sp

## sdm 1.0-89 (2020-04-22)

library(dismo)

## Loading required package: raster

library(usdm)
library(mapview)

Georreferenciación

Después cargamos las ocurrencias del oso y revisamos cuántas hay:

Oso <- dismo::gbif(
  genus="Tremarctos",
  species="ornatus", 
  download=TRUE,
  geo=TRUE,
  sp=FALSE)

## 7310 records found

## 0-300-600-900-1200-1500-1800-2100-2400-2700-3000-3300-3600-3900-4200-4500-4800-5100-5400-5700-6000-6300-6600-6900-7200-7310 records downloaded

Hay 7310 ocurrencias. Con el parámetro geo sólo descargamos los datos georreferenciados, esto es porque GBIF admite ocurrencias sin georreferenciación. Con sp establecido en falso descargamos menos información. Ahora veamos la clase en que hemos descargado los datos así como sus dimensiones.

class(Oso)

## [1] "data.frame"

dim(Oso)

## [1] 7310  157

Ah por lo anterior podemos ver que se trata de la clase data.frame de R y que sus dimensionjes son 7310 filas por 157 columnas. Ahora filtremos un poco nuestro marco de datos llamado Oso; eliminemos las filas cuya varible «longitid», lon, está ausente:

# Excluir faltantes lon
lonNAs <- which(is.na(Oso$lon))
Oso <- Oso[-lonNAs,]
# Dimensionar oso ahora
dim(Oso)

## [1] 7090  157

Normalmente también se eliman las que no tiene valor en la variable «latitud». Pero revisemos:

# Faltantes latitud
which(is.na(Oso$lat))

## integer(0)

En efecto, no hay filas sin dato para latitud. Ahora le incluimos una columna al marco de datos que represente las presencias y las ausencias. Para esto creamos una columna, Tremarctos ornatus, de unos así:

Oso$Tremarctos_ornatus <- 1

En seguida eliminamos todas las columnas quye no nos interesan y dejamos sólo las que sí:

Oso<-Oso[,c("lon","lat","Tremarctos_ornatus")]
head(Oso)

##         lon      lat Tremarctos_ornatus
## 1 -75.07116 2.751928                  1
## 2 -75.07120 2.751997                  1
## 3 -75.07094 2.752069                  1
## 4 -75.07046 2.752328                  1
## 5 -75.06934 2.750614                  1
## 6 -75.07109 2.750375                  1

Para trabajar con lso paquetes debemos convertir nuestro marco de datos a un marco de datos de puntos espaciales:

# Conversion a un "special point data frame"
sp::coordinates(Oso) <- ~lon+lat
class(Oso)

## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"

head(Oso)

##   Tremarctos_ornatus
## 1                  1
## 2                  1
## 3                  1
## 4                  1
## 5                  1
## 6                  1

Aunque no se vean, las ubicaciones siguen allí.

Variables ambientales

A continuación optenemos los datos ambientales. Para esto usamos raster::getData, con lo cual podemos obtener datos geográficos de cualquier parte del mundo. El primer parámetro es el nombre del conjunto de datos. Vamos a descargar las 19 varaibles ambientales dcon una resolución de 10 min.

bio <- raster::getData(
  name = "worldclim",
  var = "bio",
  res = 10)
bio

## class      : RasterStack 
## dimensions : 900, 2160, 1944000, 19  (nrow, ncol, ncell, nlayers)
## resolution : 0.1666667, 0.1666667  (x, y)
## extent     : -180, 180, -60, 90  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## names      :  bio1,  bio2,  bio3,  bio4,  bio5,  bio6,  bio7,  bio8,  bio9, bio10, bio11, bio12, bio13, bio14, bio15, ... 
## min values :  -269,     9,     8,    72,   -59,  -547,    53,  -251,  -450,   -97,  -488,     0,     0,     0,     0, ... 
## max values :   314,   211,    95, 22673,   489,   258,   725,   375,   364,   380,   289,  9916,  2088,   652,   261, ...

RasterStack es una superposición de rásteres. Ahora miramos la colinearidad de tales variables ambientales:

vif1 <- usdm::vifstep(bio)

## Warning in summary.lm(lm(y[, i] ~ ., data = y[-i])): essentially perfect fit:
## summary may be unreliable

vif1

## 10 variables from the 19 input variables have collinearity problem: 
##  
## bio5 bio11 bio6 bio1 bio7 bio16 bio17 bio12 bio10 bio4 
## 
## After excluding the collinear variables, the linear correlation coefficients ranges between: 
## min correlation ( bio14 ~ bio8 ):  -0.04431865 
## max correlation ( bio9 ~ bio3 ):  0.8063727 
## 
## ---------- VIFs of the remained variables -------- 
##   Variables      VIF
## 1      bio2 2.270049
## 2      bio3 4.596504
## 3      bio8 2.463748
## 4      bio9 3.342689
## 5     bio13 5.152779
## 6     bio14 3.931342
## 7     bio15 2.240163
## 8     bio18 4.217596
## 9     bio19 3.647508

Entonces es con estas 10 variables con las que vamos a trabajar porque las otras correalcionana con esas 10. Excluimos pues las varaibles con alto vif.

bio.presente <- usdm::exclude(bio,vif1)
bio.presente

## class      : RasterStack 
## dimensions : 900, 2160, 1944000, 9  (nrow, ncol, ncell, nlayers)
## resolution : 0.1666667, 0.1666667  (x, y)
## extent     : -180, 180, -60, 90  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## names      : bio2, bio3, bio8, bio9, bio13, bio14, bio15, bio18, bio19 
## min values :    9,    8, -251, -450,     0,     0,     0,     0,     0 
## max values :  211,   95,  375,  364,  2088,   652,   261,  4001,  3985

Exploremos en un mapa cómo se ven las variables ambientales que hemos dejado:

plot(bio.presente[["bio8"]]) 
# Variable ambiental bio8, "Mean Temperature of Wettest Quarter"
points(Oso,cex=0.5,pch=16) # Sobreponer los puntos georreferenciados del oso

Hay un parde datos raros, por ejemplo en Japón y en Estados Unidos. Vamos a tratar de limpiar y visualizar de mejor manera lo anterior usando . Como Oso no tiene un sistema de proyección, hay que ponerle un sistema referencial de coordenadas, en este caso se poene por defecto datum = WGS84.

# Asignacion de una projeccion
raster::projection(raster()) # ráster vacío con proyección por defecto

## [1] "+proj=longlat +datum=WGS84 +no_defs"

# asignamos tal ráster a "Oso".
sp::proj4string(Oso) <- raster::projection(raster()) 
mapview::mapview(Oso)

Para exluir los datos atípicos, que pueden ser de zoológicos o smuseos, se peude usar raster::drawExtent() y `raster::crop``:

Andes<-drawExtent()

# Limitar las ocurrencias al area de estudio
Oso<-raster::crop(Oso,Andes)
# Ver el resultado en mapas cambiando color y tamaño de lso puntos
plot(bio.presente[[1]])
points(Oso,cex=0.1,pch=16,col="red")

Elaboración de modelo para el presente

Para juntar la información georreferenciada con la ambiental usamos sdm::sdmData, que crea un objeto sdmdata que almacena la especie y sus varaibles explicativas.

Para_modelo <- sdm::sdmData(Tremarctos_ornatus~.,Oso,predictors=bio.presente)

## Loading required package: gbm

## Loaded gbm 2.1.8

## Loading required package: tree

## Loading required package: mda

## Loading required package: class

## Loaded mda 0.5-2

## Loading required package: mgcv

## Loading required package: nlme

## 
## Attaching package: 'nlme'

## The following object is masked from 'package:usdm':
## 
##     Variogram

## The following object is masked from 'package:raster':
## 
##     getData

## This is mgcv 1.8-35. For overview type 'help("mgcv-package")'.

## Loading required package: glmnet

## Loading required package: Matrix

## Loaded glmnet 4.1-2

## Loading required package: earth

## Loading required package: Formula

## Loading required package: plotmo

## Loading required package: plotrix

## Loading required package: TeachingDemos

## Loading required package: rJava

## Error: package or namespace load failed for 'rJava':
##  .onLoad failed in loadNamespace() for 'rJava', details:
##   call: fun(libname, pkgname)
##   error: JAVA_HOME cannot be determined from the Registry

## Loading required package: RSNNS

## Loading required package: Rcpp

## Loading required package: randomForest

## randomForest 4.6-14

## Type rfNews() to see new features/changes/bug fixes.

## Loading required package: rpart

## Loading required package: kernlab

## 
## Attaching package: 'kernlab'

## The following objects are masked from 'package:raster':
## 
##     buffer, rotated

## Warning in proj4string(train): CRS object has comment, which is lost in output

## Warning in proj4string(train): CRS object has comment, which is lost in output

Para_modelo

## class                                 : sdmdata 
## =========================================================== 
## number of species                     :  1 
## species names                         :  Tremarctos_ornatus 
## number of features                    :  9 
## feature names                         :  bio2, bio3, bio8, ... 
## type                                  :  Presence-Only 
## has independet test data?             :  FALSE 
## number of records                     :  2528 
## has Coordinates?                      :  TRUE

Como vemos sólo hay datos de presencia, pero debemos regrsitrar las ausencias, de modo que hay que generar seudoausencias distribuidades aleatoria e independientemente de la capa. En este causo le vamos a agregar 1000 ausencias.

Para_modelo<-sdmData(Tremarctos_ornatus~.,Oso,predictors=bio.presente,bg=list(n=1000))

## Warning in proj4string(train): CRS object has comment, which is lost in output

## Warning in proj4string(train): CRS object has comment, which is lost in output

Para_modelo

## class                                 : sdmdata 
## =========================================================== 
## number of species                     :  1 
## species names                         :  Tremarctos_ornatus 
## number of features                    :  9 
## feature names                         :  bio2, bio3, bio8, ... 
## type                                  :  Presence-Background 
## has independet test data?             :  FALSE 
## number of records                     :  3528 
## has Coordinates?                      :  TRUE

Ya con eso podemos crear nuestro modelo, pero antes hay que decidir qué algoritmo usar. Veamos qué opciones tenemos, alt = FALSE para que dé sólo los genéricos.

sdm::getmethodNames(alt=FALSE)

##  [1] "bioclim"       "bioclim.dismo" "brt"           "cart"         
##  [5] "domain.dismo"  "fda"           "gam"           "glm"          
##  [9] "glmnet"        "mahal.dismo"   "mars"          "maxent"       
## [13] "maxlike"       "mda"           "mlp"           "rbf"          
## [17] "rf"            "rpart"         "svm"

Entoces hagámoslo con los datos de presencias y ausencias. Vamos a usar un representanate de cada uno de los algoritmo que hemos visto más arriba. Vamos a correrlos dos veces n=2.

modelos <- sdm::sdm(
  Tremarctos_ornatus~.,
  Para_modelo,
  methods=c("glm","svm","rf"),
  replication=c("boot"),
  n=2)
modelos

Se carga automáticamente el paquete parallel para correr tareas en paralelo. El método de partición de os tados para tener datos de prueba y de modelaje es bootstrap. AUC, Area under the curve, sirve para medir qué tan bueno es el modelo. Para explorar más los resultados se puede usar `sdm::gui, que genera una interfaz gráfica.

sdm::gui(modelos)

A través de esa interfaz se puede evaluar los modelos por curvas ROC y AUC. Recordemos que hablamos del umbral, éste representa en el mapa ceros y unos, cada modelo recomienda un umbral según índice de desempeño.

Existe lo que se llama curvas de respuesta de la especie a cada variable.

sdm::rcurve(modelos,id=1)

Esto representa la respuesta de la especie a cada varaible, el rango en que la esepcie está según cda variable. Ahora vamos a predecir las probabilidades ocurrencia de la especie en el presente:

# Prediccion (estimación de probabilidades de ocurrencia para el presente)
presente<-predict(
  modelos,
  bio.presente,
  "prediccionesp.img", # bajar en este formato
  mean=TRUE) # da el mismo peso a las dos corridas de cada modelo

Ahora graficamos. Se muestra la probabilida de ocurrencia según los modelos.

presente
plot(presente[[c(1,2,3)]])

class : RasterBrick dimensions : 900, 2160, 1944000, 3 (nrow, ncol, ncell, nlayers) resolution : 0.1666667, 0.1666667 (x, y) extent : -180, 180, -60, 90 (xmin, xmax, ymin, ymax) crs : +proj=longlat +datum=WGS84 +no_defs source : prediccionesp.img names : sp_1.m_glm.re_boot, sp_1.m_svm.re_boot, sp_1.m_rf.re_boot min values : 6.936934e-11, -3.474973e-01, -1.635192e-15 max values : 0.9999824, 1.2331330, 1.0000000 fullname : species_Tremarctos_ornatus-method_glm-replication (Mean)_bootstrap, species_Tremarctos_ornatus-method_svm-replication (Mean)_bootstrap, species_Tremarctos_ornatus-method_rf-replication (Mean)_bootstrap

Ahora procedemos a ensamblar los modelos para obtener uno solo.

ensamble.presente <- ensemble(
  modelos,bio.presente,
  "ensamblep.img",
  setting=list(method="unweighted")
  )
#plot(ensamble.presente)

El ensamblaje predice que el osos tiene el potencial de vivir en otras regiones diferentes a los Andes. Esa sería la predicción para el presente. Pero también se peuden hacer predicciones para el futuro.

Biología computacional | Ecología

Martín Guacaneme

21/6/2021