forked from SCAR/EGABIcourse19
-
Notifications
You must be signed in to change notification settings - Fork 0
/
69_SDM_Extra.Rmd
117 lines (84 loc) · 7.25 KB
/
69_SDM_Extra.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
# Información adicional sobre modelos de distribución de especies
Esta sección contiene información adicional que no fue incorporada en el taller, pero que es relevante para los participantes.
## Predicciones parciales en el espacio geográfico
En la sección de resultados de modelos de distribución de especies incluimos información sobre gráficos de dependencia parcial. Estos gráficos muestras la contribución de cada predictor considerado en el modelo a su desempeño. Vale la pena notar que es posible evaluar estas dependencias parciales en un grilla espacial, lo cual mostrará los patrones de contribución de cada predictor en el espacio físico. Esto puede ser útil al evaluar si un predictor utilizado en el modelo se comporta de la manera esperada.
En muchas aplicaciones de modelos de distribución de especies no solemos incluir variables predictoras que nos den información sobre los procesos y condiciones que afectan directamente a la especie de nuestro interés. Esto ocurre con especial regularidad en el Océano del Sur porque es un área enorme y remote con un esfuerzo de muestreo relativamente bajo. Típicamente dependemos de información obtenida de manera remota (satélites) o producida por modelos climáticos. Las variables para las que tenemos acceso no están necesariamente directamente vinculadas al comportamiento o distribución de la especie de nuestro interés. Sino que más bien dependemos de otras variables que funcionen como indicadores.
Por ejemplo, raramente tenemos acceso a información sobre la disponibilidad de presas de nuestra especie de interés. Pero solemos tener acceso a información como clorofila deriva desde imágenes satelitales que nos dan una idea sobre la producción primaria en el área donde se encuentra nuestra especie de interés. La producción primaria nos podría servir como indicador de la disponibilidad de presas para especies que están en un nivel trófico más alto. Sin embargo, antes de incluir una variable indicadora en nuestros modelos, debemos asegurarnos que esta variable contribuya información útil sobre nuestra especie de interés.
### Ejemplo
Regresemos al ejemplo que vimos en la sección de resumen. Vamos a utilizar el mismo código que escribimos en la primera parte de esta sección.
```{r spatial_partial1, message = FALSE, results = "hide"}
library(dplyr)
library(worrms)
library(robis)
library(blueant)
library(mgcv)
library(SOmap)
library(terra)
my_cmap <- if (getRversion() >= "3.6") {
hcl.colors(21, "Geyser")
} else {
c("#008585", "#359087", "#539B8A", "#6DA590", "#85AF97", "#9BBAA0", "#AEC4AA",
"#BED0B0", "#D0DCB5", "#E5E7BC", "#FBF2C4", "#F3E3B2", "#EDD59F", "#E7C68C",
"#E3B77A", "#DEA868", "#DA9857", "#D58847", "#D1773A", "#CC6530", "#C7522B")
}
## taxonomia
x <- occurrence(datasetid = "cb16377b-56a8-4d95-802d-4eec02466773")
my_species <- "Euphausia crystallorophias"
tax <- wm_records_names(name = my_species)
my_aphia_id <- tax[[1]]$valid_AphiaID
## datos sobre presencia/ausencia
xfit <- x |>
dplyr::rename(lon = "decimalLongitude", lat = "decimalLatitude") |>
group_by(lon, lat) |>
dplyr::summarize(present = any(my_aphia_id %in% aphiaID))
## datos ambientales
## pon los datos en un directorio temporal
my_data_directory <- tempdir()
data_source <- sources_sdm("Southern Ocean marine environmental data")
## busca los datos
status <- bb_get(data_source, local_file_root = my_data_directory,
verbose = TRUE)
nc_files <- Filter(function(z) grepl("\\.nc$", z), status$files[[1]]$file)
## crea un raster apilado con las capas relevantes
env_stack <- subset(stack(nc_files), c("depth", "ice_cover_mean"))
temp <- as.data.frame(extract(env_stack, xfit[, c("lon", "lat")]))
xfit <- bind_cols(xfit, temp)
fit <- gam(present ~ s(depth) + s(ice_cover_mean) + s(lon),
family = binomial, data = xfit)
```
Este modelo `gam` de presencia/ausencia es el mismo ejemplo que vimos en la sección de Resumen, con la excepción que esta vez incluimos `lon` (longitud) como un predictor. Normalmente no utilizaríamos latitud o longitud como predictores, pero en este caso solo lo hacemos para ilustrar el ejemplo de gráfico de dependencia parcial en el espacio.
Los gráficos de dependencia parcial puede ser creados de la siguiente manera:
```{r spatial_partial2}
plot(fit, pages = 1, shade = TRUE, scale = 0)
```
Podemos ver que el término parcial para longitud (panel abajo a la izquierda) tiene ondulaciones, en especial alrededor de `-60` donde luce algo raro. Sin embarho es difícil interpretar este gráfico porque sabemos que esta variable no nos proporciona información que se relevante desde un punto de vista ecológico.
Hagamos un mapa del término parcial para ver cómo varía en el espacio.
```{r spatial_partial3}
## creamos una grilla vacia
xpred <- expand.grid(
lon = seq(from = floor(min(xfit$lon)),
to = ceiling(max(xfit$lon)), by = 0.25),
lat = seq(from = floor(min(xfit$lat)),
to = ceiling(max(xfit$lat)), by = 0.25))
## incluimos datos ambientales en la grilla que acabamos de crear
xpred <- bind_cols(as.data.frame(xpred),
as.data.frame(extract(env_stack,
xpred[, c("lon", "lat")])))
## termimnos de prediccion parcial
sp_pred <- predict(fit, newdata = xpred, type = "terms")
```
Usamos `predict(..., type = "terms")` para obtener la contribución de cada término en el modelo cuando predicimos la presencia de especies basada en la información en la grilla `xpred`. La opción `type = "terms"` funciona solamente con unos pocos modelos, como GLMs y GAMs, pero existen métodos similares que pueden ser aplicados a otro tipo de modelos.
Ahora podemos hacer un gráfico de la contribución de `lon` en el área que estamos modelando:
```{r spatial_partial4}
xpred <- cbind(xpred, sp_pred)
pr <- rast(xpred[, c("lon", "lat", "s(lon)")], type = "xyz")
## el termino "s(lon)" introduce un spline de suavizado (smooth) a la variable
# entre parentesis
## grafica el resultado
projection(pr) <- "+proj=longlat +datum=WGS84"
p <- SOmap_auto(x = pr, bathy = pr)
p$bathy_palette <- my_cmap
p
```
En este gráfico podemos ver que el término `lon` simplemente permite que el modelo prediga probabilidad más alta de presencia en ciertas bandas de longitud. Este término no da información que sea relevante desde un punto de vista ecológico, así que probablemente debamos descartar este término del modelo. Sería mejor que encontremos otra variable indicadora que nos de información más relevante sobre la ecología de la especie que estamos intentando modelar.
Algo que debemos considerar es que el término `lon` parece resultar en mejores predicciones, pero debemos considerar lo que significa desde un punto de vista ecológico. Es posible que este término sea de gran importancia para mejorar el modelo, pero esto puede ocurrir de casualidad. No porque explique la ecología o comportamiento de nuestra especie de interés, sino porque explique algún otro patrón en nuestro set de datos de entrenamiento. Por esta razón debemos estar atentos de incluir variables que expliquen procesos ecológicos de otra manera corremos el riesgo de tener un modelo que nos de resultados sin relevancia ecológica.