¿Cómo crear un modelo de distribución consensuado?

Para crear un modelo de distribución consensuado, puedes combinar múltiples modelos de distribución de especies individuales (SDMs) en un solo modelo. Aquí hay los pasos generales para hacerlo:

Antes de comenzar: creación del espacio de trabajo

Antes de empezar, es necesario instalar los paquetes de R requeridos. Por favor, realiza una sola ejecución de instalación por computadora. Recuerda que la instalación se realiza con la función install.packages() y el nombre del paquete entre comillas. Por ejemplo, install.packages(“maptools”). Luego, puedes cargar los paquetes en tu sesión de R utilizando la función library(), como se muestra en el siguiente ejemplo: library(maptools). A continuación los paquetes:

library(dismo)
library(maptools)
library(sp)
library(tidyverse)
library(raster)
library(rgdal)
library(sf)
library(randomForest)
library(kernlab)
library(bioclim)

nota: Quiza, deberas actualizar Java para el paquete “randomForest”.

Elige un espacio de trabajo

Para trabajar en R, es importante elegir un espacio de trabajo o una carpeta. Se recomienda usar la función setwd para establecer el directorio de trabajo, pero hay otras formas de hacerlo. El siguiente es un ejemplo de cómo utilizar setwd:

setwd("D:/biato/wikiaves/clase_conservacion")
###### Este ultimo vamos a usarlo para otra cosa pero... debe ser el mismo de arriba 
first_folder<-"D:/biato/wikiaves/clase_conservacion"

nota: Aquí hay un ejemplo. Es recomendable que utilices el mismo nombre final de carpeta, “clase_conservacion”, para almacenar toda tu información.

Importar puntos de datos

Debe guardar su data frame como archivo “.csv” y con el tipo CSV UTF (separado por comas). A continuación, encontrará un ejemplo con una planta de la familia Lauraceae, pero se recomienda que le de un nombre diferente a su archivo de data frame; por favor, asegúrese de utilizar nombres simples y sin espacios o caracteres especiales.

#nota que no use espacios ni asentos en el nombre de mi .csv
obs.data <- read.csv(file = "D:/biato/wikiaves/clase_conservacion/Comino_Crespo.csv", header = T, sep = ",", dec = ".")

Luego usa la funcion “head(obs.data, 4)” para ver si estan de forma correcta. Debes ver algo como lo siguiente,

Ejemplo de como deben estar tus datos.
gbifID datasetKey occurrenceID kingdom phylum class order family genus species infraspecificEpithet taxonRank scientificName verbatimScientificName verbatimScientificNameAuthorship countryCode locality stateProvince occurrenceStatus individualCount publishingOrgKey decimalLatitude decimalLongitude coordinateUncertaintyInMeters coordinatePrecision elevation elevationAccuracy depth depthAccuracy eventDate day month year taxonKey speciesKey basisOfRecord institutionCode collectionCode catalogNumber recordNumber identifiedBy dateIdentified license rightsHolder recordedBy typeStatus establishmentMeans lastInterpreted mediaType issue
931017814 23243164-c132-444d-a3a4-34662abb8305 UdeA:HUA:97603 Plantae Tracheophyta Magnoliopsida Laurales Lauraceae Aniba Aniba perutilis NA SPECIES Aniba perutilis Hemsl. Aniba perutilis Hemsl. CO Hacienda Peñaflor. cerro Caliza. quebrada Caño Seco Antioquia PRESENT 1 cccff716-2694-4209-9f9e-2f7a484465a0 6.272225 -74.68528 1295 NA 300 0 NA NA 1995-12-13T00:00:00 13 12 1995 4180651 4180651 PRESERVED_SPECIMEN Universidad de Antioquia (UdeA) HUA 97603 52 Sánchez C. 2011-01-01T00:00:00 CC0_1_0 Gíl N. NA NA 2023-01-24T15:50:55.178Z COORDINATE_ROUNDED;INSTITUTION_MATCH_FUZZY;COLLECTION_MATCH_FUZZY
931009904 23243164-c132-444d-a3a4-34662abb8305 UdeA:HUA:6270 Plantae Tracheophyta Magnoliopsida Laurales Lauraceae Aniba Aniba perutilis NA SPECIES Aniba perutilis Hemsl. Aniba perutilis Hemsl. CO Corregimiento Providencia. valle del río Anorí. entre Dos Bocas y Anorí. cerca a Planta Providencia. 26 km S and 23 km W air of Zaragoza Antioquia PRESENT 1 cccff716-2694-4209-9f9e-2f7a484465a0 7.216667 -75.05056 30214 NA 550 150 NA NA 1977-08-01T00:00:00 1 8 1977 4180651 4180651 PRESERVED_SPECIMEN Universidad de Antioquia (UdeA) HUA 6270 828 Callejas R. 1992-01-01T00:00:00 CC0_1_0 Shepherd J. NA NA 2023-01-24T15:51:09.700Z COORDINATE_ROUNDED;INSTITUTION_MATCH_FUZZY;COLLECTION_MATCH_FUZZY
931009853 23243164-c132-444d-a3a4-34662abb8305 UdeA:HUA:6249 Plantae Tracheophyta Magnoliopsida Laurales Lauraceae Aniba Aniba perutilis NA SPECIES Aniba perutilis Hemsl. Aniba perutilis Hemsl. CO Corregimiento Providencia. valle del río Anorí. entre Dos Bocas y Anorí. cerca a Planta Providencia. 26 km S and 23 km W air of Zaragoza Antioquia PRESENT 1 cccff716-2694-4209-9f9e-2f7a484465a0 7.216667 -75.05056 30214 NA 550 150 NA NA 1977-08-01T00:00:00 1 8 1977 4180651 4180651 PRESERVED_SPECIMEN Universidad de Antioquia (UdeA) HUA 6249 821 Sánchez C. 2010-08-01T00:00:00 CC0_1_0 Shepherd J. NA NA 2023-01-24T15:51:09.648Z COORDINATE_ROUNDED;INSTITUTION_MATCH_FUZZY;COLLECTION_MATCH_FUZZY
931003789 23243164-c132-444d-a3a4-34662abb8305 UdeA:HUA:32317 Plantae Tracheophyta Magnoliopsida Laurales Lauraceae Aniba Aniba perutilis NA SPECIES Aniba perutilis Hemsl. Aniba perutilis Hemsl. CO Vereda Berlín Antioquia PRESENT 1 cccff716-2694-4209-9f9e-2f7a484465a0 5.966667 -74.83333 1486 NA 275 25 NA NA 1986-11-23T00:00:00 23 11 1986 4180651 4180651 PRESERVED_SPECIMEN Universidad de Antioquia (UdeA) HUA 32317 1183 Hoyos S. CC0_1_0 Hoyos S.;Alzate N. NA NA 2023-01-24T15:51:05.795Z COORDINATE_ROUNDED;INSTITUTION_MATCH_FUZZY;COLLECTION_MATCH_FUZZY

nota: Como puedes ver, el data-frame de “GBif” contiene mucha información que, aunque es significativa, no es relevante para nuestro análisis. Por lo tanto, debemos limpiar los datos y utilizar solo las columnas necesarias.

Limpiar y procesar nuestro data-frame

Mejor debemos limpiar y eliminar alguna información de nuestro data frame para evitar errores y problemas en el código. Solo vamos a tomar las columnas con información geográfica; debido a un análisis de una sola especie en nuestro ejemplo. Se supone que todas las coordenadas son de una especie. Para esto, seleccionaremos las columnas “decimalLongitude” a “decimalLatitude”, que son las que tienen la longitud y latitud, respectivamente. Tambien, limpiaremos los NA y cambiaremos el nombre de las columnas por unas mas faciles (Long y Lat, respectivamente).

acgeo2 <- obs.data %>%  subset( !is.na(decimalLongitude ) & !is.na(decimalLatitude)) %>% 
  transform(as.numeric(decimalLongitude)) %>% transform(as.numeric(decimalLatitude)) %>% 
  select(decimalLongitude, decimalLatitude) %>% 
  rename(Long=decimalLongitude) %>% 
  rename(Lat=decimalLatitude)

######### tu solo corre esto, el se encargara de limpiar los datos y colocar nombres mas simples

tus datos se veran de la siguiente manera…

Ejemplo de tus datos luego de la limpieza.
Long Lat
-74.68528 6.272225
-75.05056 7.216667
-75.05056 7.216667
-74.83333 5.966667

nota: Asegúrate de revisar tus datos utilizando la función “head(acgeo2, 4)”.

Descarga las capas climáticas de la base de datos WWF y el raster de composición política de Colombia

Descargaremos algunos raster Biogeográficos de internet para nuestro ejercicio. Hay muchas formas de utilizar más raster en nuestro análisis y proceso, pero debido a ser un ejemplo simple, solo usaremos esos.

files <- list.files(path=paste(system.file(package="dismo"),
                               '/ex', sep=''), pattern='grd', full.names=T)

Bioma <- stack(files)

nota: Nuestro objeto “Bioma” tiene todos los raster importados que vamos a utilizar.

Descarga de la división política del país

El paquete ‘raster’ incluye la división política de algunos países. En el caso colombiano, tiene el país (capa 0), departamentos (capa 1) y regiones (capa 2).

Colombia0 <- getData('GADM' , country="COL", level=0)
Colombia1 <- getData('GADM' , country="COL", level=1)

plot(Colombia0, main="Administrativas Colombia nivel 0")

plot(Colombia1, main="Adm. Colombia nivel 1")

plot(Colombia0, xlim=c(-90,-60), ylim=c(-5,10), axes=TRUE, col="light yellow")
points(acgeo2$Lon, acgeo2$Lat, col='orange', pch=20, cex=0.75)
points(acgeo2$Lon, acgeo2$Lat, col='red', cex=0.75)

**nota*: Preguntate, ¿estan todos los puntos dentro de Colombia? En caso que no, los descargaste mal.

Continúa procesando

Necesitamos ajustar el raster bio climático a una extensión geográfica específica, en este caso Colombia. Para ello, crearemos una zona de trabajo geográfica.

e <- extent(-81.841530, -66.87033, -4.228429 , 15.91247) #Colombia extencion

Y ahora lo “cortamos”…

Bioma <- crop(x = Bioma, y = e) ##this funtion cuts "Bioma" rasters by "e" extent
predictors2 <- stack(Bioma)

plot(predictors2)

Quieres ver los datos por cada punto, bien has lo siguiente…

presvals <- raster::extract(predictors2, acgeo2) 

head(presvals, 3) 
##      bio1 bio12 bio16 bio17 bio5 bio6 bio7 bio8 biome
## [1,]  252  2664   919   334  308  198  110  249     1
## [2,]  232  3921  1364   365  290  171  119  229     1
## [3,]  232  3921  1364   365  290  171  119  229     1

Problema frecuente

En algunos casos dependiendo la version del computador o el R, etc… Se pueden presentar problemas en el acceso de las coordenas antes suministradas como coordenadas y no numeros. Por tanto realizamos lo siguiente, aclarando a R que son coordenadas.

Long<-as.numeric(as.matrix(acgeo2[1]))
Lat<-as.numeric(as.matrix(acgeo2[2]))
coord2 <- cbind(Long, Lat)
coord2<-as.data.frame(coord2)
acgeo3<-raster::coordinates(coord2)

presvals <-  raster::extract(predictors2, acgeo3)

Creando puntos de ausencia (adsence points)

Existen muchas maneras de modelar la distribución de una especie, y la mayoría de ellas pueden ser creadas por una mezcla de información de presencia / ausencia. Para lograr nuestro objetivo, debemos crear algunos puntos de ausencia pseudo-aleatorios.

set.seed(0)
backgr <- randomPoints(predictors2, 500)
absvals <-  raster::extract(predictors2, backgr)
pb <- c(rep(1, nrow(presvals)), rep(0, nrow(absvals)))
sdmdata <- data.frame(cbind(pb, rbind(presvals, absvals)))
sdmdata[,'biome'] <- as.factor(sdmdata[,'biome'])

pred_nf <- dropLayer(predictors2, 'biome')

set.seed(0)
group <- kfold(acgeo2, 5)
pres_train <- acgeo3[group != 1, ]
pres_test <- acgeo3[group == 1, ]

####         
set.seed(100)
backg <- randomPoints(pred_nf, n=10000, ext=e, extf = 1.25)

Estos pueden apreciarse de la siguiente manera:

colnames(backg) = c('Lon', 'Lat')
group <- kfold(backg, 5)
backg_train <- backg[group != 1, ]
backg_test <- backg[group == 1, ]

r <- raster(pred_nf, 1)
plot(!is.na(r), col=c('white', 'light grey'), legend=FALSE)
plot(e, add=TRUE, col='red', lwd=2)
points(backg_train, pch='-', cex=0.5, col='yellow')
points(backg_test, pch='-',  cex=0.5, col='black')
points(pres_train, pch= '+', col='green')
points(pres_test, pch='+', col='blue')

Modelar la distribución de especies

La distribución de especies es fundamental en la conservación biológica, ya que nos permite entender cómo las especies interactúan con su entorno y cómo pueden verse afectadas por cambios ambientales. Conocer la distribución de una especie también nos ayuda a identificar áreas clave para su conservación y tomar medidas para protegerlas. Además, modelar la distribución de especies nos permite predecir cómo una especie puede reaccionar a futuros cambios climáticos y en su hábitat, lo cual es esencial para planificar la conservación a largo plazo.

A continuación, mostraremos cómo realizar algunos modelos de distribución. Al finalizar el proceso, crearemos un consenso de los modelos y exportaremos el resultado.

Método clásico Bioclim (climate-envelope-model)

El modelo de sobre-clima es un tipo de modelo de distribución de especies que predice la distribución de una especie basado en sus requisitos climáticos. El modelo funciona creando una “envoltura climática” que define el rango de condiciones climáticas adecuadas para la especie, basadas en su distribución observada. La envoltura climática se representa típicamente por una combinación de variables bioclimáticas, como la temperatura, la precipitación y la altitud, que se utilizan para definir la tolerancia de la especie a diferentes condiciones ambientales.

bc <- bioclim(pred_nf, pres_train)
#### numero de variables climaticas que afectan la distribucion estan en verde
plot(bc, a=1, b=2, p=0.85)

ef <- evaluate(pres_test, backg_test, bc, pred_nf)
ef
## class          : ModelEvaluation 
## n presences    : 90 
## n absences     : 157 
## AUC            : 0.9803255 
## cor            : 0.8807402 
## max TPR+TNR at : 0.04133646
tr <- threshold(ef, 'spec_sens')
tr
## [1] 0.04133646
pb <- predict(pred_nf, bc, ext=e, progress='')
pb
## class      : RasterLayer 
## dimensions : 40, 30, 1200  (nrow, ncol, ncell)
## resolution : 0.5, 0.5  (x, y)
## extent     : -82, -67, -4, 16  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : layer 
## values     : 0, 0.8425414  (min, max)

Grafica:

par(mfrow=c(1,2))
plot(pb, main='Bioclim, raw values')
plot(Colombia0, add=TRUE, border='dark grey')
plot(pb > tr, main='presence/absence')
plot(Colombia0, add=TRUE, border='dark grey')
points(pres_train, pch='+')

Método de dominio

En el contexto del modelado de distribución, el método de dominio se refiere a un enfoque específico para definir la extensión geográfica, o dominio, del modelo. El método de dominio se utiliza para limitar el área de estudio a la gama relevante de una especie y para excluir áreas donde la especie es poco probable que ocurra. Esto ayuda a mejorar la precisión y eficiencia del modelo de distribución, evitando la extrapolación más allá del rango de la especie y centrándose en los recursos computacionales en las áreas más importantes.

dm <- domain(pred_nf, pres_train)
ef <- evaluate(pres_test, backg_test, dm, pred_nf)
ef
## class          : ModelEvaluation 
## n presences    : 90 
## n absences     : 157 
## AUC            : 0.9559448 
## cor            : 0.712134 
## max TPR+TNR at : 0.9393305

Grafica:

##Grafica
pd = predict(pred_nf, dm, ext=e, progress='')
par(mfrow=c(1,2))
plot(pd, main='Domain, raw values')
plot(Colombia0, add=TRUE, border='dark grey')
tr <- threshold(ef, 'spec_sens')
plot(pd > tr, main='presence/absence')
plot(Colombia0, add=TRUE, border='dark grey')
points(pres_train, pch='*')

La distancia de Mahalanobis

La distancia de Mahalanobis es una medida de la distancia entre un punto y una distribución. Es comúnmente utilizada en análisis estadístico para evaluar la diferencia entre una observación y un grupo de observaciones, teniendo en cuenta la covarianza entre las variables.

Se define como la distancia entre un punto y una distribución medida en desviaciones estándar, utilizando la matriz de covarianza de la distribución para ajustar por la correlación entre variables. A menudo se utiliza en reconocimiento de patrones y aprendizaje automático, especialmente en la clasificación de datos multivariados.

mm <- mahal(pred_nf, pres_train)


ef <- evaluate(pres_test, backg_test, mm, pred_nf)
ef
## class          : ModelEvaluation 
## n presences    : 90 
## n absences     : 157 
## AUC            : 0.9872611 
## cor            : 0.2968733 
## max TPR+TNR at : 0.9999
pm = predict(pred_nf, mm, ext=e, progress='')
par(mfrow=c(1,2))
pm[pm < -10] <- -10

Grafica:

##Graphic
plot(pm, main='Mahalanobis distance')
plot(Colombia0, add=TRUE, border='dark grey')

tr <- threshold(ef, 'spec_sens')
plot(pm > tr, main='presence/absence')
plot(Colombia0, add=TRUE, border='dark grey')
points(pres_train, pch='+')

Modelado por regresion glm

data3 <- backg_train 
colnames(data3) <- colnames(pres_train)
rbind(pres_train, data3)
backg_train<-data3 
train <- rbind(pres_train, backg_train)
pb_train <- c(rep(1, nrow(pres_train)), rep(0, nrow(backg_train)))
envtrain <-  raster::extract(predictors2, train)
envtrain <- data.frame( cbind(pa=pb_train, envtrain) )
envtrain[,'biome'] = factor(envtrain[,'biome'], levels=1:14)
head(envtrain)
##   pa bio1 bio12 bio16 bio17 bio5 bio6 bio7 bio8 biome
## 1  1  252  2664   919   334  308  198  110  249     1
## 2  1  232  3921  1364   365  290  171  119  229     1
## 3  1  252  2664   919   334  308  198  110  249     1
## 4  1  264  2498   922   372  324  211  113  260     1
## 5  1  252  2664   919   334  308  198  110  249     1
## 6  1  262  3743  1323   345  321  206  115  258     1
testpres <- data.frame(  raster::extract(predictors2, pres_test) )
testbackg <- data.frame(  raster::extract(predictors2, backg_test) )
testpres[ ,'biome'] = factor(testpres[ ,'biome'], levels=1:14)
testbackg[ ,'biome'] = factor(testbackg[ ,'biome'], levels=1:14)



##Generalized Linear Models
gm1 <- glm(pa ~ bio1 + bio5 + bio6 + bio7 + bio8 + bio12 + bio16 + bio17,
           family = binomial(link = "logit"), data=envtrain)


coef(gm1)
##  (Intercept)         bio1         bio5         bio6         bio7         bio8 
##  9.025910183 -0.024883294  2.067536524 -1.974100826 -2.062222389 -0.087604751 
##        bio12        bio16        bio17 
##  0.007446127 -0.019399999 -0.011623238
gm2 <- glm(pa ~ bio1+bio5 + bio6 + bio7 + bio8 + bio12 + bio16 + bio17,
           family = gaussian(link = "identity"), data=envtrain)
evaluate(testpres, testbackg, gm1)
## class          : ModelEvaluation 
## n presences    : 90 
## n absences     : 157 
## AUC            : 0.9504246 
## cor            : 0.760687 
## max TPR+TNR at : 1.274852
ge2 <- evaluate(testpres, testbackg, gm2)
ge2
## class          : ModelEvaluation 
## n presences    : 90 
## n absences     : 157 
## AUC            : 0.9526893 
## cor            : 0.7768101 
## max TPR+TNR at : 0.7905402
pg <- predict(predictors2, gm2, ext=e)
par(mfrow=c(1,2))

Grafica:

##Graphic
plot(pg, main='GLM/gaussian, raw values')
plot(Colombia0, add=TRUE, border='dark grey')

tr <- threshold(ge2, 'spec_sens')
plot(pg > tr, main='presence/absence')
plot(Colombia0, add=TRUE, border='dark grey')
points(pres_train, pch='+')
points(backg_train, pch='-', cex=0.25)

Modelado por regresion glm

Los métodos de aprendizaje automático son un tipo de inteligencia artificial que implican entrenar algoritmos para hacer predicciones basadas en datos de entrada. Las técnicas de aprendizaje automático, como árboles de decisión, bosques aleatorios y máquinas de vectores de soporte, pueden utilizarse para modelar las distribuciones de especies relacionando la presencia o ausencia de una especie con variables ambientales.

Por otro lado, Maxent es un tipo de algoritmo de aprendizaje automático diseñado específicamente para modelar la distribución de especies. Maxent utiliza un enfoque de entropía máxima para modelar la relación entre la presencia de una especie y las variables ambientales, teniendo en cuenta la incertidumbre en las predicciones del modelo.

Importante Deberas intalas un paquete extra con la actualizacion de Java.

#install.packages("rJava") #Corre la instalacion solo una vez OJOOOOOO
library(rJava)

Ahora si continuemos!!

xm <- maxent(predictors2, pres_train, factors='biome')
plot(xm)

ef <- evaluate(pres_test, backg_test, xm, predictors2)
ef
## class          : ModelEvaluation 
## n presences    : 90 
## n absences     : 157 
## AUC            : 0.9398797 
## cor            : 0.7830463 
## max TPR+TNR at : 0.8959265
px <- predict(predictors2, xm, ext=e, progress='')
par(mfrow=c(1,2))

Grafica:

## Grafica
plot(px, main='Maxent, raw values')
plot(Colombia0, add=TRUE, border='dark grey')

tr <- threshold(ef, 'spec_sens')
plot(px > tr, main='presence/absence')
plot(Colombia0, add=TRUE, border='dark grey')
points(pres_train, pch='+')

Random Forest

Random Forest es una técnica de aprendizaje automático utilizada para el modelado predictivo. Construye un conjunto de árboles de decisión y combina sus resultados para hacer una predicción final para nuevas entradas. Además, ofrece una evaluación de la importancia de las características y está equipado para procesar datos faltantes y desequilibrados.

model <- pa ~ bio1 + bio5 + bio6 + bio7 + bio8 + bio12 + bio16 + bio17
rf1 <- randomForest(model, data=envtrain)
model <- factor(pa) ~ bio1 + bio5 + bio6 + bio7 + bio8 + bio12 + bio16 + bio17
rf2 <- randomForest(model, data=envtrain)
rf3 <- randomForest(envtrain[,1:8], factor(pb_train))
erf <- evaluate(testpres, testbackg, rf1)
erf
## class          : ModelEvaluation 
## n presences    : 90 
## n absences     : 157 
## AUC            : 0.9950814 
## cor            : 0.9613489 
## max TPR+TNR at : 0.5890039
pr <- predict(predictors2, rf1, ext=e)

Grafica:

par(mfrow=c(1,2))
plot(pr, main='Random Forest, regression')
plot(Colombia0, add=TRUE, border='dark grey')
tr <- threshold(erf, 'spec_sens')
plot(pr > tr, main='presence/absence')
plot(Colombia0, add=TRUE, border='dark grey')
points(pres_train, pch='+')
points(backg_train, pch='-', cex=0.25)

Máquinas de Vectores de Soporte (Support Vector Machines;SVM)

El método de Máquinas de Vectores de Soporte (SVM) también puede ser utilizado en modelado de distribución de especies. ksvm en R permite la aplicación de SVM en problemas de modelado de distribución de especies, utilizando técnicas de kernel para resolver problemas no lineales y obtener una mejor precisión en la predicción de la distribución de especies en un área dada.

svm <- ksvm(pa ~ bio1+bio5+bio6+bio7+bio8+bio12+bio16+bio17, data=envtrain)
esv <- evaluate(testpres, testbackg, svm)
esv
## class          : ModelEvaluation 
## n presences    : 90 
## n absences     : 157 
## AUC            : 0.9881458 
## cor            : 0.8938387 
## max TPR+TNR at : 0.1064286
ps <- predict(predictors2, svm, ext=e)

Grafica:

### Graphic
par(mfrow=c(1,2))
plot(ps, main='Support Vector Machine')
plot(Colombia0, add=TRUE, border='dark grey')
tr <- threshold(esv, 'spec_sens')
plot(ps > tr, main='presence/absence')
plot(Colombia0, add=TRUE, border='dark grey')
points(pres_train, pch='+')
points(backg_train, pch='-', cex=0.25)

Consenso de motodos

El método de consenso consiste en combinar los resultados de varios modelos diferentes para obtener una predicción más precisa y robusta. Se puede crear un método de consenso mediante la combinación de las predicciones de varios modelos individuales utilizando técnicas como la media, la votación o la combinación ponderada. Es importante seleccionar cuidadosamente los modelos individuales que se van a combinar y evaluar su rendimiento antes de crear el método de consenso.

models <- stack(pb, pd, pm, pg, pr, ps )
names(models) <- c("bioclim", "domain", "mahal", "glm", "Random Forest", "SVM")

Grafica de comparacion:

### Graphic
plot(models)

Tarea: Que diferencias puedes encontrar? Meditalo y analiza criticamente cual es la diferencia entre lo que modelamos, la realida y cuales son las mayores interrogantes?

Ahora si podemos crear el consenso final

m <- mean(models)
m<- mask(m,Colombia0 )

Grafica final:

plot(m, main='average score')
points(pres_train, pch='.')
points(backg_train, pch='-', cex=0.25)
plot(Colombia0, add=TRUE, border='dark grey')

Exportamos resultados

Ahora podemos exportar los resultados para un análisis más exhaustivo de la distribución. Analicemos si los puntos de presencia coinciden con la distribución prevista.

auc <- sapply(list(ge2, erf, esv), function(x) x@auc)
w <- (auc-0.5)^2
plot(m, main='weighted mean of three models')
points(pres_train, pch='+')
points(backg_train, pch='-', cex=0.25)

Vamos a crear una carperta llamada distribucion…

file <- "distribucion"

dir.create(file.path(first_folder, file )) 
## Warning in dir.create(file.path(first_folder, file)):
## 'D:\biato\wikiaves\clase_conservacion\distribucion' already exists

Ve a tu carpeta “clase_conservacion”, en ella ya podras apreciar una carpeta llamada “distribucion”. en ella vamos a guardar nuestro raster con la informacion de distribucion.

setwd("D:/biato/wikiaves/clase_conservacion/distribucion")

#### por favor ten cuidado con los nombres de las carpetasssssssssss, en tu pc deben ser distintas

pol <- rasterToPolygons(m>0.4,function(x) x == 1,dissolve=T) ##depends on your specie you would ##increase or decrease the porcentage 
raster::shapefile(pol,
                  "Specie_fin.shp", 
                  overwrite=TRUE ) ### Ojo con los nombres de las carpetas, en

Visualizacion de resultado final

Como paso extra vamos a visualizar la distribucion en R y jugar un poco con esta. Para el siguiente paso debemos intalar y ejecutar mas paquetes. Recuerda que la instalación se realiza con la función install.packages() y el nombre del paquete entre comillas. Por ejemplo, install.packages(“maptools”). Luego, puedes cargar los paquetes en tu sesión de R utilizando la función library(), como se muestra en el siguiente ejemplo: library(maptools). A continuación los paquetes:

library(htmlwidgets)
library(webshot)
library(leaflet)
library(mapview)

##########Lo siguiente es para fijar nuestra carpeta como de grabado de outputs
setwd("D:/biato/wikiaves/clase_conservacion/distribucion")
getwd()
## [1] "D:/biato/wikiaves/clase_conservacion/distribucion"

A continuacion crearemos una serie de condiciones con el fin de crear un mapa interactivo en R, desde el cual podremos interactuar de forma mas natural con nuestra distribucion.

##### Recuerda lo de las carpetas, aqui debes colocar el orden correcto y nombres


wqp_sites <- st_read("D:/biato/wikiaves/clase_conservacion/distribucion")
## Reading layer `Specie_fin' from data source 
##   `D:\biato\wikiaves\clase_conservacion\distribucion' using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 1 field
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -76.5 ymin: 1.5 xmax: -73 ymax: 7.5
## Geodetic CRS:  GCS_unknown
nhd_wms_url <- "https://server.arcgisonline.com/ArcGIS/rest/services/World_Physical_Map/MapServer/tile/{z}/{y}/{x}"
#######
##### Aca creamos una funcion para algunas especificaciones en nuestro mapa interactivo
extractCoords <- function(sp.df)
{
  results <- list()
  for(i in 1:length(sp.df@polygons[[1]]@Polygons))
  {
    results[[i]] <- sp.df@polygons[[1]]@Polygons[[i]]@coords
  }
  results <- Reduce(rbind, results)
  results
}

Ahora usaremos un paquete llamado “leaflet” el cual usa una interfase html para disponer mapas de forma didactica.

## OGR data source with driver: ESRI Shapefile 
## Source: "D:\biato\wikiaves\clase_conservacion\distribucion\Specie_fin.shp", layer: "Specie_fin"
## with 1 features
## It has 1 fields