Modelamiento macroecológico
Velasco-Vinasco J.A.
2024-11-04
En esta sección haremos un ejercicio de modelación de riqueza de especies o modelo macroecológico (MEM en breve) en función de un conjunto de variables bioclimáticas y evaluaremos la capacidad predictiva del modelo con datos por fuera de la muestra y haremos algunas proyecciones a futuro usando escenarios de cambio climático.
Descargue los datos de aquí: https://www.dropbox.com/sh/9ks083qx3f0jm5m/AADOj-uvf3tUEGtHaBQ_1Kosa?dl=0
##rgdal, rgeos and maptools replaced by sf, stars and terra
library(raster)## Loading required package: splibrary(sp)
library(sf)## Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1; sf_use_s2() is TRUElibrary(stars)## Loading required package: abindlibrary(terra)## terra 1.7.83#library(maptools)
library(rgdal)## Please note that rgdal will be retired during October 2023,
## plan transition to sf/stars/terra functions using GDAL and PROJ
## at your earliest convenience.
## See https://r-spatial.org/r/2023/05/15/evolution4.html and https://github.com/r-spatial/evolution
## rgdal: version: 1.6-7, (SVN revision 1203)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.0.4, released 2020/01/28
## Path to GDAL shared files: /usr/share/gdal
## GDAL binary built with GEOS: TRUE
## Loaded PROJ runtime: Rel. 6.3.1, February 10th, 2020, [PJ_VERSION: 631]
## Path to PROJ shared files: /usr/share/proj
## Linking to sp version:2.1-4
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading sp or rgdal.##
## Attaching package: 'rgdal'## The following object is masked from 'package:terra':
##
## projectlibrary(dismo)
library(XML)
library(maps)
library(letsR)
library(plyr)##
## Attaching package: 'plyr'## The following object is masked from 'package:maps':
##
## ozonelibrary(groupdata2)
#library(rgeos)
library(classInt)
library(Metrics)
library(RColorBrewer)
library(rcompanion)##
## Attaching package: 'rcompanion'## The following object is masked from 'package:Metrics':
##
## accuracylibrary(plotKML)## plotKML version 0.8-3 (2022-06-05)## URL: https://github.com/Envirometrix/plotKML/library(lmtest)## Loading required package: zoo##
## Attaching package: 'zoo'## The following object is masked from 'package:terra':
##
## time<-## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericlibrary(spdep)## Loading required package: spData## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`library(lme4)## Loading required package: Matrix##
## Attaching package: 'lme4'## The following object is masked from 'package:raster':
##
## getDatalibrary(reshape2)
library(usdm)
setwd("/home2/taller_R/")Cargar shapefile con areas de distribución (extent of occurrence EOO) de todas las especies de ranas de la familia Centrolenidae (ranas de cristal).
crs.geo <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84") # geographical, datum WGS84
mapas <-st_read("./Shapefiles/Centrolenidae_IUCN.shp")## Reading layer `Centrolenidae_IUCN' from data source
## `/home2/taller_R/Shapefiles/Centrolenidae_IUCN.shp' using driver `ESRI Shapefile'
## Simple feature collection with 158 features and 28 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -100.3641 ymin: -29.32361 xmax: -37.46495 ymax: 19.56507
## Geodetic CRS: GCS_unknownmapas## Simple feature collection with 158 features and 28 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -100.3641 ymin: -29.32361 xmax: -37.46495 ymax: 19.56507
## Geodetic CRS: GCS_unknown
## First 10 features:
## id_no binomial presence origin seasonal compiler
## 1 54917 Centrolene hesperium 1 1 1 Kelsey Neam
## 2 54959 Cochranella euhystrix 1 1 1 Kelsey Neam
## 3 54979 Nymphargus phenax 1 1 1 Kelsey Neam
## 4 54979 Nymphargus phenax 1 1 1 Steve Best
## 5 54931 Centrolene peristictum 1 1 1 IUCN
## 6 78457419 Centrolene sabini 1 1 1 Kelsey Neam
## 7 54933 Centrolene pipilatum 1 1 1 IUCN
## 8 78458868 Cochranella guayasamini 1 1 1 Kelsey Neam
## 9 78903218 Chimerella corleone 1 1 1 Kelsey Neam
## 10 78926514 Hyalinobatrachium anachoretus 1 1 1 Kelsey Neam
## year_
## 1 2018
## 2 2017
## 3 2018
## 4 2018
## 5 2008
## 6 2017
## 7 2008
## 8 2017
## 9 2018
## 10 2018
## citation
## 1 IUCN (International Union for Conservation of Nature), Conservation International & NatureServe
## 2 IUCN (International Union for Conservation of Nature), Conservation International & NatureServe
## 3 IUCN (International Union for Conservation of Nature), Conservation International & NatureServe
## 4 IUCN (International Union for Conservation of Nature), Conservation International & NatureServe
## 5 IUCN (International Union for Conservation of Nature), Conservation International & NatureServe.
## 6 IUCN (International Union for Conservation of Nature), Conservation International & NatureServe.
## 7 IUCN (International Union for Conservation of Nature), Conservation International & NatureServe.
## 8 IUCN (International Union for Conservation of Nature), Conservation International & NatureServe
## 9 IUCN (International Union for Conservation of Nature), Conservation International & NatureServe
## 10 IUCN (International Union for Conservation of Nature), Conservation International & NatureServe
## source dist_comm island subspecies subpop legend tax_comm
## 1 <NA> <NA> <NA> <NA> <NA> Extant (resident) <NA>
## 2 <NA> <NA> <NA> <NA> <NA> Extant (resident) <NA>
## 3 <NA> <NA> <NA> <NA> <NA> Extant (resident) <NA>
## 4 <NA> <NA> <NA> <NA> <NA> Extant (resident) <NA>
## 5 <NA> <NA> <NA> <NA> <NA> Extant (resident) <NA>
## 6 <NA> <NA> <NA> <NA> <NA> Extant (resident) <NA>
## 7 <NA> <NA> <NA> <NA> <NA> Extant (resident) <NA>
## 8 <NA> <NA> <NA> <NA> <NA> Extant (resident) <NA>
## 9 <NA> <NA> <NA> <NA> <NA> Extant (resident) <NA>
## 10 <NA> <NA> <NA> <NA> <NA> Extant (resident) <NA>
## kingdom phylum class order_ family genus category
## 1 ANIMALIA CHORDATA AMPHIBIA ANURA CENTROLENIDAE Centrolene EN
## 2 ANIMALIA CHORDATA AMPHIBIA ANURA CENTROLENIDAE Cochranella CR
## 3 ANIMALIA CHORDATA AMPHIBIA ANURA CENTROLENIDAE Nymphargus EN
## 4 ANIMALIA CHORDATA AMPHIBIA ANURA CENTROLENIDAE Nymphargus EN
## 5 ANIMALIA CHORDATA AMPHIBIA ANURA CENTROLENIDAE Centrolene VU
## 6 ANIMALIA CHORDATA AMPHIBIA ANURA CENTROLENIDAE Centrolene VU
## 7 ANIMALIA CHORDATA AMPHIBIA ANURA CENTROLENIDAE Centrolene EN
## 8 ANIMALIA CHORDATA AMPHIBIA ANURA CENTROLENIDAE Cochranella LC
## 9 ANIMALIA CHORDATA AMPHIBIA ANURA CENTROLENIDAE Chimerella DD
## 10 ANIMALIA CHORDATA AMPHIBIA ANURA CENTROLENIDAE Hyalinobatrachium EN
## marine terrestial freshwater SHAPE_Leng SHAPE_Area
## 1 f t t 0.25286444 0.0025440442
## 2 f t t 0.12643222 0.0012720221
## 3 f t t 0.09823846 0.0007679566
## 4 f t t 0.09823846 0.0007679566
## 5 f t t 16.27430080 0.9418965598
## 6 f t t 0.10123419 0.0008155085
## 7 f t t 0.73770207 0.0117393643
## 8 f t t 2.62224604 0.2668484008
## 9 f t t 0.10662779 0.0009047232
## 10 f t t 0.21783713 0.0018953110
## species geometry
## 1 Centrolene_hesperium MULTIPOLYGON (((-79.08572 -...
## 2 Cochranella_euhystrix MULTIPOLYGON (((-79.06595 -...
## 3 Nymphargus_phenax MULTIPOLYGON (((-73.59147 -...
## 4 Nymphargus_phenax MULTIPOLYGON (((-73.91681 -...
## 5 Centrolene_peristictum MULTIPOLYGON (((-76.19387 4...
## 6 Centrolene_sabini MULTIPOLYGON (((-71.6086 -1...
## 7 Centrolene_pipilatum MULTIPOLYGON (((-77.89176 -...
## 8 Cochranella_guayasamini MULTIPOLYGON (((-76.0394 -7...
## 9 Chimerella_corleone MULTIPOLYGON (((-76.2782 -6...
## 10 Hyalinobatrachium_anachoretus MULTIPOLYGON (((-77.64369 -...plot(mapas)## Warning: plotting the first 9 out of 28 attributes; use max.plot = 28 to plot
## all
st_transform(mapas,crs = crs.geo)## Simple feature collection with 158 features and 28 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -100.3641 ymin: -29.32361 xmax: -37.46495 ymax: 19.56507
## Geodetic CRS: GEOGCRS["unknown",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]],
## ID["EPSG",6326]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8901]],
## CS[ellipsoidal,2],
## AXIS["longitude",east,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433,
## ID["EPSG",9122]]],
## AXIS["latitude",north,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433,
## ID["EPSG",9122]]]]
## First 10 features:
## id_no binomial presence origin seasonal compiler
## 1 54917 Centrolene hesperium 1 1 1 Kelsey Neam
## 2 54959 Cochranella euhystrix 1 1 1 Kelsey Neam
## 3 54979 Nymphargus phenax 1 1 1 Kelsey Neam
## 4 54979 Nymphargus phenax 1 1 1 Steve Best
## 5 54931 Centrolene peristictum 1 1 1 IUCN
## 6 78457419 Centrolene sabini 1 1 1 Kelsey Neam
## 7 54933 Centrolene pipilatum 1 1 1 IUCN
## 8 78458868 Cochranella guayasamini 1 1 1 Kelsey Neam
## 9 78903218 Chimerella corleone 1 1 1 Kelsey Neam
## 10 78926514 Hyalinobatrachium anachoretus 1 1 1 Kelsey Neam
## year_
## 1 2018
## 2 2017
## 3 2018
## 4 2018
## 5 2008
## 6 2017
## 7 2008
## 8 2017
## 9 2018
## 10 2018
## citation
## 1 IUCN (International Union for Conservation of Nature), Conservation International & NatureServe
## 2 IUCN (International Union for Conservation of Nature), Conservation International & NatureServe
## 3 IUCN (International Union for Conservation of Nature), Conservation International & NatureServe
## 4 IUCN (International Union for Conservation of Nature), Conservation International & NatureServe
## 5 IUCN (International Union for Conservation of Nature), Conservation International & NatureServe.
## 6 IUCN (International Union for Conservation of Nature), Conservation International & NatureServe.
## 7 IUCN (International Union for Conservation of Nature), Conservation International & NatureServe.
## 8 IUCN (International Union for Conservation of Nature), Conservation International & NatureServe
## 9 IUCN (International Union for Conservation of Nature), Conservation International & NatureServe
## 10 IUCN (International Union for Conservation of Nature), Conservation International & NatureServe
## source dist_comm island subspecies subpop legend tax_comm
## 1 <NA> <NA> <NA> <NA> <NA> Extant (resident) <NA>
## 2 <NA> <NA> <NA> <NA> <NA> Extant (resident) <NA>
## 3 <NA> <NA> <NA> <NA> <NA> Extant (resident) <NA>
## 4 <NA> <NA> <NA> <NA> <NA> Extant (resident) <NA>
## 5 <NA> <NA> <NA> <NA> <NA> Extant (resident) <NA>
## 6 <NA> <NA> <NA> <NA> <NA> Extant (resident) <NA>
## 7 <NA> <NA> <NA> <NA> <NA> Extant (resident) <NA>
## 8 <NA> <NA> <NA> <NA> <NA> Extant (resident) <NA>
## 9 <NA> <NA> <NA> <NA> <NA> Extant (resident) <NA>
## 10 <NA> <NA> <NA> <NA> <NA> Extant (resident) <NA>
## kingdom phylum class order_ family genus category
## 1 ANIMALIA CHORDATA AMPHIBIA ANURA CENTROLENIDAE Centrolene EN
## 2 ANIMALIA CHORDATA AMPHIBIA ANURA CENTROLENIDAE Cochranella CR
## 3 ANIMALIA CHORDATA AMPHIBIA ANURA CENTROLENIDAE Nymphargus EN
## 4 ANIMALIA CHORDATA AMPHIBIA ANURA CENTROLENIDAE Nymphargus EN
## 5 ANIMALIA CHORDATA AMPHIBIA ANURA CENTROLENIDAE Centrolene VU
## 6 ANIMALIA CHORDATA AMPHIBIA ANURA CENTROLENIDAE Centrolene VU
## 7 ANIMALIA CHORDATA AMPHIBIA ANURA CENTROLENIDAE Centrolene EN
## 8 ANIMALIA CHORDATA AMPHIBIA ANURA CENTROLENIDAE Cochranella LC
## 9 ANIMALIA CHORDATA AMPHIBIA ANURA CENTROLENIDAE Chimerella DD
## 10 ANIMALIA CHORDATA AMPHIBIA ANURA CENTROLENIDAE Hyalinobatrachium EN
## marine terrestial freshwater SHAPE_Leng SHAPE_Area
## 1 f t t 0.25286444 0.0025440442
## 2 f t t 0.12643222 0.0012720221
## 3 f t t 0.09823846 0.0007679566
## 4 f t t 0.09823846 0.0007679566
## 5 f t t 16.27430080 0.9418965598
## 6 f t t 0.10123419 0.0008155085
## 7 f t t 0.73770207 0.0117393643
## 8 f t t 2.62224604 0.2668484008
## 9 f t t 0.10662779 0.0009047232
## 10 f t t 0.21783713 0.0018953110
## species geometry
## 1 Centrolene_hesperium MULTIPOLYGON (((-79.08572 -...
## 2 Cochranella_euhystrix MULTIPOLYGON (((-79.06595 -...
## 3 Nymphargus_phenax MULTIPOLYGON (((-73.59147 -...
## 4 Nymphargus_phenax MULTIPOLYGON (((-73.91681 -...
## 5 Centrolene_peristictum MULTIPOLYGON (((-76.19387 4...
## 6 Centrolene_sabini MULTIPOLYGON (((-71.6086 -1...
## 7 Centrolene_pipilatum MULTIPOLYGON (((-77.89176 -...
## 8 Cochranella_guayasamini MULTIPOLYGON (((-76.0394 -7...
## 9 Chimerella_corleone MULTIPOLYGON (((-76.2782 -6...
## 10 Hyalinobatrachium_anachoretus MULTIPOLYGON (((-77.64369 -...st_crs(mapas)## Coordinate Reference System:
## User input: GCS_unknown
## wkt:
## GEOGCRS["GCS_unknown",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]],
## ID["EPSG",6326]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["Degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["longitude",east,
## ORDER[1],
## ANGLEUNIT["Degree",0.0174532925199433]],
## AXIS["latitude",north,
## ORDER[2],
## ANGLEUNIT["Degree",0.0174532925199433]]]library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:plyr':
##
## arrange, count, desc, failwith, id, mutate, rename, summarise,
## summarize## The following objects are masked from 'package:terra':
##
## intersect, union## The following objects are masked from 'package:raster':
##
## intersect, select, union## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionmapas2<-mapas %>% group_by(binomial) %>% summarise()Generar un mapa de riqueza y una matriz de presencia-ausencia (PAM) usando los mapas individuales.
PAM_centro <-lets.presab(mapas2, xmn=-100.4, xmx=-37.5, ymn=-29.3, ymx=19.6, resol=1,
remove.cells=T, remove.sp=T, show.matrix=FALSE,
crs=CRS("+proj=longlat +ellps=WGS84 +datum=WGS84"), cover=0, presence=NULL,
origin=NULL, seasonal=NULL, count=FALSE)plot(PAM_centro,main="Riqueza de especies")
PAM_centro##
## Class: PresenceAbsence
## Number of species: 148
## Number of cells: 514
## Resolution: 1, 1 (x, y)str(PAM_centro)## List of 3
## $ Presence_and_Absence_Matrix: num [1:514, 1:150] -96.9 -95.9 -94.9 -96.9 -95.9 -94.9 -93.9 -92.9 -91.9 -99.9 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : NULL
## .. ..$ : chr [1:150] "Longitude(x)" "Latitude(y)" "Celsiella revocata" "Celsiella vozmedianoi" ...
## $ Richness_Raster :S4 class 'SpatRaster' [package "terra"]
## $ Species_name : chr [1:148] "Celsiella revocata" "Celsiella vozmedianoi" "Centrolene acanthidiocephalum" "Centrolene altitudinale" ...
## - attr(*, "class")= chr "PresenceAbsence"Extraer raster
rr <- PAM_centro$Richness_RasterCargar un mapa de países (shapefile)
am <-readOGR("./Shapefiles/america_completo.shp")## Warning in readOGR("./Shapefiles/america_completo.shp"): OGR support is
## provided by the sf and terra packages among others## Warning in ogrListLayers(dsn = dsn): OGR support is provided by the sf and
## terra packages among others## Warning in ogrInfo(dsn = dsn, layer = layer, encoding = encoding, use_iconv =
## use_iconv, : OGR support is provided by the sf and terra packages among others## Warning in ogrFIDs(dsn = dsn, layer = layer): OGR support is provided by the sf
## and terra packages among others## Warning in OGRSpatialRef(dsn, layer, morphFromESRI = morphFromESRI, dumpSRS =
## dumpSRS, : OGR support is provided by the sf and terra packages among others## Warning in ogrListLayers(dsn): OGR support is provided by the sf and terra
## packages among others## Warning in ogrFIDs(dsn = dsn, layer = layer): OGR support is provided by the sf
## and terra packages among others## OGR data source with driver: ESRI Shapefile
## Source: "/home2/taller_R/Shapefiles/america_completo.shp", layer: "america_completo"
## with 5763 features
## It has 7 fieldsprojection(am)<-crs.geocolor<-hcl.colors(n=100,rev = T,palette="viridis")
color<-c("white",color)
plot(rr, main="Riqueza de especies",col=color)
plot(am, add=T)
Subir capas bioclimáticas actuales y futuras
bios.act <- stack(list.files(path="./Bioclimaticos/actual/", pattern = "tif", full.names = T))
bios.fut <- stack(list.files(path="./Bioclimaticos/rcp85_2070/bc85bi70/", pattern = "tif", full.names = T))
projection(bios.act)<-crs.geo
projection(bios.fut)<-crs.geoExtraer datos de variables ambientales para cada celda con datos de riqueza
PAM_centro_bios <- lets.addvar(PAM_centro, bios.act, fun=mean, onlyvar = TRUE)
PAM_centro_biosF <- lets.addvar(PAM_centro, bios.fut, fun=mean, onlyvar = TRUE)Preparar los conjuntos de datos de riqueza de especies y variables bioclimáticas
PAM_DF <- as.data.frame(PAM_centro$Presence_and_Absence_Matrix)
a <- rowSums(PAM_DF[,c(-1:-2)])
DF_SR <- PAM_DF[,c(1:2)]
DF_SR["Riqueza"] <- a
DF_SR_Bios <- round(cbind(DF_SR, PAM_centro_bios),1)
DF_SR_BiosFut <- round(cbind(DF_SR, PAM_centro_biosF),1)Construir el modelo de regresión. Esta parte es muy importante porque es necesario tener una buena calibración y que el modelo tenga buena capacidad de predicción de la riqueza por fuera de la muestra. Primero vamos a generar una partición aleatoria de los datos. Hay otro tipo de particiones estructuradas que se pueden usar (e.g., stratified, block, checkerboard, véase el siguiente enlace: http://cran.nexr.com/web/packages/ENMeval/vignettes/ENMeval-vignette.html#partition).
partitions <- partition(DF_SR_Bios, p=0.7)
train <- partitions[[1]]
test <- partitions[[2]]
trainDF <- as.data.frame(train)
testDF <- as.data.frame(test)Ajustar el modelo. Usaremos una regresión lineal con el método de mínimos cuadrados ordinarios (OLS en inglés). Este método hace una estimación de los parámetros minimizando la suma de residuales al cuadrado. Es el método más sencillo.
m.ols <- lm(Riqueza~., data=trainDF[,c(-1,-2)])
summary(m.ols)##
## Call:
## lm(formula = Riqueza ~ ., data = trainDF[, c(-1, -2)])
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.9783 -1.4688 -0.1985 1.1361 12.5116
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -27.941536 10.248717 -2.726 0.006737 **
## bio_01_mean -0.194774 0.207610 -0.938 0.348824
## bio_02_mean -0.484414 0.127040 -3.813 0.000163 ***
## bio_03_mean 0.627975 0.139426 4.504 9.19e-06 ***
## bio_04_mean -0.016301 0.008534 -1.910 0.056947 .
## bio_05_mean -0.576730 2.824328 -0.204 0.838319
## bio_06_mean 0.816060 2.823255 0.289 0.772720
## bio_07_mean 1.043996 2.819702 0.370 0.711427
## bio_08_mean -0.039085 0.064273 -0.608 0.543522
## bio_09_mean 0.079438 0.059622 1.332 0.183636
## bio_10_mean 0.478318 0.290410 1.647 0.100476
## bio_11_mean -0.629217 0.394685 -1.594 0.111817
## bio_12_mean 0.003765 0.001708 2.204 0.028179 *
## bio_13_mean 0.066231 0.013043 5.078 6.31e-07 ***
## bio_14_mean 0.013822 0.035888 0.385 0.700363
## bio_15_mean 0.005197 0.017303 0.300 0.764070
## bio_16_mean -0.036208 0.005173 -7.000 1.38e-11 ***
## bio_17_mean -0.026636 0.013277 -2.006 0.045628 *
## bio_18_mean 0.017046 0.002283 7.466 7.03e-13 ***
## bio_19_mean 0.005097 0.001507 3.381 0.000806 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.582 on 339 degrees of freedom
## Multiple R-squared: 0.5728, Adjusted R-squared: 0.5489
## F-statistic: 23.93 on 19 and 339 DF, p-value: < 2.2e-16Generamos la predicción usando el conjunto de datos que se dejaron por fuera (fuera de la muestra).
pred.ols <- predict(m.ols, testDF)
# Comparación entre observaciones y predicciones.
par(mfrow=c(1,1))
hist(testDF$Riqueza, main="Observaciones")
hist(pred.ols, main="Predicciones (fuera de la muestra)")
nsim <- 100
mfit_ols <- list() # guarda cada modelo ajustado
rsq_ols <- list() # guarda cada R2 calculado
rmse_ols <- list() # guarda la raíz del error cuadrático medio RMSE
pred_ols <- list() # valores predichos fuera de la muestra
for (i in 1:nsim){
#
partitions <- partition(DF_SR_Bios, p=0.7)
train <- partitions[[1]]
test <- partitions[[2]]
trainDF <- as.data.frame(train)
testDF <- as.data.frame(test)
# model fit
m.ols <- lm(Riqueza~., data=trainDF[,c(-1:-2)])
mfit_ols[[i]] <- m.ols
# predictions
pred.ols <- predict(m.ols, testDF[,c(-1:-2)])
pred_ols[[i]] <- pred.ols
# metrics
rss <- sum((pred.ols - testDF$Riqueza) ^ 2) ## suma de cuadrados residuales
tss <- sum((testDF$Riqueza - mean(testDF$Riqueza)) ^ 2) ## suma total de cuadrados
rsq_ols[[i]] <- 1 - rss/tss # R2
rmse_ols[[i]] <- rmse(testDF$Riqueza, pred.ols) # raíz del error cuadrático medio RMSE
}Hacemos un histograma de los 100 valores de R2 y RMSE
bb <- cbind(as.numeric(rsq_ols), as.numeric(rmse_ols))
colnames(bb) <- c("R2", "RMSE")
v.ols <- as.data.frame(bb)
par(mfrow=c(1,2))
hist(v.ols$R2, main="R2")
hist(v.ols$RMSE, main="RMSE")
Ahora vamos a ajustar un modelo linear generalizado. Es más flexible que
el OLS y permite acomodar mejorar la variación de los datos porque la
variable de respuesta puede tener un tipo de distribución diferente a la
Gaussiana (no normal).
# GLM
nsim <- 100
# to run models across a sample of 100 spatial folds
mfit_glm <- list() # model fitted
rsq_glm <- list() # R squared
rmse_glm <- list() # Root Mean Squared Error
pred_glm <- list() # predicciones fuera de la muestra
for (i in 1:nsim){
#
partitions <- partition(DF_SR_Bios, p=0.7)
train <- partitions[[1]]
test <- partitions[[2]]
trainDF <- as.data.frame(train)
testDF <- as.data.frame(test)
# model fit
m.glm <- glm(Riqueza~., data = trainDF[,c(-1:-2)], family=poisson())
mfit_glm[[i]] <- m.glm
# predictions
pred.glm <- predict(m.glm, newdata=testDF[,c(-1:-2)], type="response", se.fit=FALSE)
pred_glm[[i]] <- pred.glm
# metrics
rss <- sum((pred.glm - testDF$Riqueza) ^ 2) ## suma de cuadrados residuales
tss <- sum((testDF$Riqueza - mean(testDF$Riqueza)) ^ 2) ## suma total de cuadrados
rsq_glm[[i]] <- 1 - rss/tss # R2
rmse_glm[[i]] <- rmse(testDF$Riqueza, pred.glm) # raíz del error cuadrático medio RMSE
}Hacemos el histograma de los 100 valores de R2 y RMSE y comparamos el GLM con el OLS
bb <- cbind(as.numeric(rsq_glm), as.numeric(rmse_glm))
colnames(bb) <- c("R2", "RMSE")
v.glm <- as.data.frame(bb)
par(mfrow=c(1,2))
hist(v.ols$R2, main="R2 OLS")
hist(v.ols$RMSE, main="RMSE OLS")
par(mfrow=c(1,2))
hist(v.glm$R2, main="R2 GLM")
hist(v.glm$RMSE, main="RMSE GLM")
Vamos a generar unos escenarios futuros de riqueza de especies. Recuerde las diferencias entre pronóstico, predicción y proyección. Cuando hablamos de cambio climático el término adecuado es proyección o escenario futuro.
# Proyectar a futuro
options(scipen=999) # quitamos la notación científica
m.ols <- lm(Riqueza~., data=DF_SR_Bios[,c(-1:-2)])
m.glm <- glm(Riqueza~., data=DF_SR_Bios[,c(-1:-2)], family=poisson())
pred.ols.act <- round(predict(m.ols, newdata=DF_SR_Bios[,c(-1:-2)], type="response", se.fit=FALSE),1)
pred.glm.act <- round(predict(m.glm, newdata=DF_SR_Bios[,c(-1:-2)], type="response", se.fit=FALSE),1)
pred.ols.fut <- round(predict(m.ols, newdata=DF_SR_BiosFut[,c(-1:-2)], type="response", se.fit=FALSE),1)
pred.glm.fut <- round(predict(m.glm, newdata=DF_SR_BiosFut[,c(-1:-2)], type="response", se.fit=FALSE),1)Comparamos las ganancias y pérdidas potenciales de especies a futuro usando el modelo OLS y el GLM
pp.ols <- pred.ols.fut - pred.ols.act # negativo = pérdidas, positivo = ganancias
pp.glm <- pred.glm.fut - pred.glm.act # negativo = pérdidas, positivo = gananciasHacemos un dataframe con las predicciones futuras
DF_SR_Bios["pred_ols_act"] <- pred.ols.act
DF_SR_Bios["pred_ols_fut"] <- pred.ols.fut
DF_SR_Bios["pred_glm_act"] <- pred.glm.act
DF_SR_Bios["pred_glm_fut"] <- pred.glm.futConvertimos las predicción en formato raster
DF_SR_Bios_sp <- DF_SR_Bios
DF_SR_Bios_sp <- plyr::rename(DF_SR_Bios_sp, replace=c("Longitude(x)"="x", "Latitude(y)"="y"))
coordinates(DF_SR_Bios_sp)<- ~ x + y
class(DF_SR_Bios_sp)## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"obs.r <- vect2rast(DF_SR_Bios_sp, fname = "Riqueza", cell.size=1)
obs.r <- raster(obs.r)
projection(obs.r)<-crs.geo
ols.act.r <- vect2rast(DF_SR_Bios_sp, fname = "pred_ols_act", cell.size=1)
ols.act.r <- raster(ols.act.r)
projection(ols.act.r)<-crs.geo
ols.fut.r <- vect2rast(DF_SR_Bios_sp, fname = "pred_ols_fut", cell.size=1)
ols.fut.r <- raster(ols.fut.r)
projection(ols.fut.r)<-crs.geo
glm.act.r <- vect2rast(DF_SR_Bios_sp, fname = "pred_glm_act", cell.size=1)
glm.act.r <- raster(glm.act.r)
projection(glm.act.r)<-crs.geo
glm.fut.r <- vect2rast(DF_SR_Bios_sp, fname = "pred_glm_fut", cell.size=1)
glm.fut.r <- raster(glm.fut.r)
projection(glm.fut.r)<-crs.geoCalculamos las diferencias entre lo actual y lo futuro
dif.ols <- ols.fut.r - ols.act.r # negativo = pérdidas, positivo = ganancias
plot(dif.ols,col=color)
plot(am, add=T)
dif.glm <- glm.fut.r - glm.act.r # negativo = pérdidas, positivo = ganancias
plot(dif.glm,col=color)
plot(am,add=T)
perd.riq.ols <- reclassify(dif.ols, c(0, Inf, 0))
perd.riq.glm <- reclassify(dif.glm, c(0, Inf, 0))
gan.riq.ols <- reclassify(dif.ols, c(-Inf, 0, 0))
gan.riq.glm <- reclassify(dif.glm, c(-Inf, 0, 0))Hacemos unos mapas bonitos
library(viridis)## Loading required package: viridisLite##
## Attaching package: 'viridis'## The following object is masked from 'package:maps':
##
## unemppar(mfrow=c(2,3))
plot(ols.act.r, main="OLS predicción", col=viridis(30))
plot(am, add=T, lwd=0.5)
plot(perd.riq.ols, main="OLS pérdidas", col=viridis(30))
plot(am, add=T, lwd=0.5)
plot(gan.riq.ols, main="OLS ganancias", col=viridis(30))
plot(am, add=T, lwd=0.5)
plot(glm.act.r, main="GLM predicción", col=viridis(30))
plot(am, add=T, lwd=0.5)
plot(perd.riq.glm, main="GLM pérdidas", col=viridis(30))
plot(am, add=T, lwd=0.5)
plot(gan.riq.glm, main="GLM ganancias", col=viridis(30))
plot(am, add=T, lwd=0.5)
Estos modelos son muy básicos porque solo incorporan variables
bioclimáticas y la riqueza de especies está influenciada por otros
factores. Aquí una literatura adicional recomendada al respecto en la
que hemos trabajado con varios colegas para entender que determina la
riqueza de especies en diferentes escalas y con diferentes grupos
taxonómicos.
García-Rodríguez, A., Martínez, P.A., Oliveira, B.F., Velasco, J.A., Pyron, A & Costa, G.C. 2021. Amphibian speciation rates support a general role of mountains as biodiversity pumps. The American Naturalist. 198: E68-E79. https://doi.org/10.1086/715500
García-Rodriguez, A., Velasco, J.A, Villalobos, F & Parra, G. 2021. Effects of evolutionary time, speciation rates and local abiotic conditions on the origin and maintenance of amphibian montane diversity. Global Ecology and Biogeography 30: 674-684. https://doi.org/10.1111/geb.13249
Ochoa-Ochoa, L.,Mejía-Domínguez, N.R., Velasco, J.A., Dimitrov, D & Marske, K.A. 2020. Dimensions of Amphibian alpha diversity in the New World. Journal of Biogeography 47: 2293-2302.
Ochoa-Ochoa, L. M., Mejía-Domínguez, N.R.,Velasco, J.A, Marskea, K.A. & Rahbek, C. 2019. Amphibian functional diversity is related to precipitation stability in the New World. Global Ecology and Biogeography. 28: 1219-1229.
Velasco, J.A., Villalobos, F., Diniz-Filho, JAF., Algar, A.C., Flores-Villela, O., Kohler, G., Poe, S., Martínez-Meyer, E. 2018. Climatic and evolutionary factors shaping geographical gradients of species richness in Anolis lizards. Biological Journal of the Linnean Society 123: 615-627.
Eso es todo por ahora!