¿Quién soy?

Chris es Profesor en el Instituto de Estudios Marinos y Antárticos de la Universidad de Tasmania. Chris y su equipo en el laboratorio Seascapemodels trabajan en la conservación de ecosistemas oceánicos y la gestión sostenible de pesquerías. Su equipo utiliza avances en enfoques de modelado estadístico para sintetizar datos ecológicos e informar la toma de decisiones ambientales.

La carrera de Chris se construyó con las herramientas de código abierto de R. Él devuelve ese favor creando recursos gratuitos de R y enseñando talleres en todas partes.

Introducción

R es el lenguaje de programación líder para el modelado ecológico por buenas razones. Ser gratuito y de código abierto ciertamente ayuda. Tener fortalezas en visualización de datos también ayuda. Y debido a estas características, R ahora tiene un enorme ecosistema de paquetes contribuidos por usuarios que permiten aplicaciones de modelado sofisticadas.

Este ecosistema de paquetes de R es creado por una enorme comunidad de usuarios de R, muchos de los cuales son líderes en su campo desarrollando modelos estadísticos de vanguardia y herramientas de ciencia de datos. Esto significa que si conoces R, puedes acceder a herramientas de vanguardia y combinarlas de nuevas maneras.

Aunque hay otros lenguajes que sobresalen en computación científica, R domina el mercado en modelado estadístico ecológico, el tema de este curso.

Hasta hace poco, la mayoría de los usuarios de R preparaban sus datos fuera de R (por ejemplo, en Excel o Arc GIS) y luego los leían en R para el modelado. Pero R ahora también tiene paquetes de SIG y mapeo eficientes y fáciles de usar. Esto significa que puedes crear todo tu flujo de trabajo de análisis, desde la descarga de datos hasta la visualización, en R.

Pero comenzar un proyecto en R puede ser desalentador para nuevos usuarios. Hay una curva de aprendizaje pronunciada para la codificación, y hay tantas opciones de paquetes que es fácil tener parálisis de decisión.

Así que en este curso vamos a proporcionar una introducción a algunos enfoques comunes de modelado en R. También aprenderemos cómo construir un flujo de trabajo eficiente. Nuestros objetivos hoy son:

Revisar la capacidad de R para análisis de datos, gráficos y mapeo
Aprender a construir flujos de trabajo eficientes y repetibles para modelado predictivo
Aprender a ejecutar algunos modelos
Aprender a visualizar datos espaciales y resultados de modelos

Este curso es adecuado para personas que tienen alguna experiencia con R. No es un curso para principiantes, pero los usuarios con un poco de experiencia en R pueden seguir las primeras secciones.

Métodos que cubriremos

En este curso revisaremos:

Mapeo en R con archivos de forma y rásteres (usando los paquetes modernos sf y terra)
Modelos lineales generalizados
Modelos aditivos generalizados

Lo que necesitarás

Asegúrate de tener una versión reciente de R (disponible en el sitio web de CRAN) y Rstudio instalado (la versión gratuita de escritorio está bien).

Necesitarás instalar estos paquetes:

install.packages(c("tmap", "tidyverse",
                   "sf", "corrplot",
                   "patchwork", "visreg"))

También necesitarás descargar los datos para el curso.

Caso de estudio: Pez loro de cabeza de joroba, ‘Topa’ en las Islas Salomón

Los peces loro de cabeza de joroba (Bolbometopon muricatum) son una especie enigmática de peces tropicales. Los adultos de esta especie se caracterizan por una gran protuberancia en su frente que los machos usan para exhibirse y luchar durante la reproducción. La determinación del sexo para esta especie es desconocida, pero es probable que un individuo tenga el potencial de desarrollarse como macho o hembra al alcanzar la madurez.

Los adultos viajan en escuelas y consumen algas mordiendo trozos de coral y en el proceso literalmente defecan arena limpia. Debido a su gran tamaño, hábito de formar cardúmenes y madurez tardía, son susceptibles a la sobrepesca, y muchas poblaciones están en declive.

Su ciclo de vida se caracteriza por la migración desde arrecifes lagunares como juveniles (ver imagen a continuación) hasta hábitats de plataforma de arrecife y arrecifes expuestos como adultos. Los juveniles en etapas tempranas son carnívoros y se alimentan de zooplancton, y luego se transforman en herbívoros a una edad temprana.

Imagen: Ciclo de vida del pez loro de cabeza de joroba. Imagen de E. Stump y obtenida de Hamilton et al. 2017.

Hasta mediados de la década de 2010, el hábitat para las postlarvas y juveniles en asentamiento era un misterio. Sin embargo, el patrón de migración desde la costa hacia mar adentro a lo largo de su ciclo de vida conocido sugiere que las primeras etapas bentónicas (‘reclutas’) podrían ocurrir en hábitats de arrecifes cercanos a la costa.

Los hábitats de arrecifes cercanos a la costa son susceptibles a la degradación por mala calidad del agua, lo que genera preocupaciones de que esta especie también pueda estar en declive debido a la contaminación. Pero la falta de datos de las primeras etapas de vida dificulta una mayor exploración de este problema.

En este curso analizaremos la primera encuesta que reveló las preferencias de hábitat de las primeras etapas juveniles del pez loro de cabeza de joroba. Estos datos fueron analizados por Hamilton et al. 2017 y Brown and Hamilton 2018.

En la década de 2010, Rick Hamilton (The Nature Conservancy) dirigió una serie de encuestas en los hábitats de arrecifes cercanos a la costa de la provincia de Kia, Islas Salomón. El objetivo era buscar el hábitat de reclutamiento para los peces loro juveniles de cabeza de joroba. Estas encuestas fueron motivadas por la preocupación de las comunidades locales en Kia de que los topa (el nombre local para los peces loro de cabeza de joroba) están en declive.

En las encuestas, los buzos nadaron transectos estandarizados y buscaron peces loro juveniles de cabeza de joroba en hábitats cercanos a la costa, a menudo a lo largo del borde de los manglares. En total, encuestaron 49 sitios en toda Kia.

Estas encuestas fueron aún más desafiantes debido a la presencia de cocodrilos en el hábitat de manglares de la región. Por lo tanto, estos datos son increíblemente valiosos.

La tala en la región de Kia ha causado problemas de calidad del agua que pueden afectar los hábitats coralinos cercanos a la costa. Durante la tala, los troncos se transportan desde la tierra a las barcazas en ‘estanques de troncos’. Un estanque de troncos es un área de manglares que se arrasa con bulldozers para permitir la transferencia de troncos a las barcazas. Como puedes imaginar, los estanques de troncos son muy fangosos. Este daño crea una escorrentía significativa de sedimentos que puede asfixiar y matar los hábitats coralinos.

Rick y el equipo encuestaron arrecifes cerca de estanques de troncos y en áreas que no tenían tala. Solo encontraron reclutas de peces loro de cabeza de joroba escondiéndose en especies de coral ramificado.

En este curso, primero preguntaremos si la presencia de reclutas de peces loro de cabeza de joroba está relacionada con la cobertura de especies de coral ramificado. Luego desarrollaremos un modelo estadístico para analizar la relación entre la contaminación de los estanques de troncos y los reclutas de peces loro de cabeza de joroba, y usaremos este modelo para predecir los impactos de la contaminación en los peces loro de cabeza de joroba en toda la región de Kia.

Los datos y el código para los análisis originales están disponibles en mi sitio de github. En este curso utilizaremos versiones simplificadas de los datos originales. Agradecemos a Rick Hamilton por proporcionar los datos para este curso.

Planificando nuestro proyecto

Una parte importante de la programación en R es tener un flujo de trabajo organizado. Estar organizado es importante en cualquier cosa más compleja que un programa de R de un solo paso. Organizar tu flujo de trabajo requiere planificación anticipada y adherirse a algunas estrategias. En este curso seguiremos un flujo de trabajo simple.

Hay múltiples beneficios de un flujo de trabajo organizado. Tu código será más transparente y repetible, lo que beneficiará tanto a tu futuro yo como a otros investigadores. También significa que es menos probable que cometas errores y el código es más fácil de depurar cuando los haces. Finalmente, es posible que quieras compartir el código públicamente para que otras personas puedan repetir tu trabajo.

Primero, querrás identificar tu pregunta de investigación. Si tu pregunta es clara, entonces puedes ceñirte a lo que necesitas hacer, y no terminar yendo por caminos que no llevan a ninguna parte creando código que no usarás. Una vez que tengas esa pregunta, recomiendo considerar cinco pasos en tu flujo de trabajo:

1. Recopilar los datos

Necesitarás obtener datos, a menudo de repositorios de datos en línea. Incluso si has recopilado observaciones de especies tú mismo, es probable que necesites obtener covariables ambientales para la predicción de fuentes como oficinas meteorológicas, repositorios oceanográficos o modelos climáticos.

2. Manipulación de datos

Los datos deben ser leídos en R y ‘manipulados’ en la estructura de datos apropiada para los paquetes de modelado que utilizarás. Solo una advertencia, algunos paquetes requieren diferentes estructuras de datos, ¡así que asegúrate de saber lo que quieres! Para un modelo espacial, este paso puede involucrar mucho SIG (complejo). Afortunadamente, R tiene buenos paquetes para SIG.

3. Visualización de datos

Antes de comenzar, siempre uso las poderosas herramientas de visualización de datos de R para explorar los datos. Esto te da una comprensión más profunda de lo que modelarás, puede ayudar a evitar fallas conceptuales. Además, es posible que desees guardar algunos gráficos de datos para tu publicación. En este curso usaremos R markdown para hacer un informe sobre nuestros datos que se pueda compartir fácilmente con colaboradores.

4. Modelado

Finalmente, llegamos al modelado. Aquí es donde pasaremos la mayor parte del tiempo en este curso.

5. Visualización de modelos

Una forma poderosa de comunicar tus modelos es haciendo más visualizaciones de datos. En este curso usaremos visualizaciones de datos para ver qué dicen los modelos sobre los impulsores ambientales importantes de la abundancia de peces.

Una nota final, lo anterior puede ser un proceso iterativo. Pero organiza tus scripts de R como se indicó anteriormente y será mucho más fácil.

Comencemos

Eso es todo el trasfondo, comencemos. Trabajaremos a través de cada paso del flujo de trabajo anterior.

Configurando tu carpeta de proyecto

Inicia Rstudio. Recomiendo configurar cada proyecto como un proyecto de Rstudio. Haz clic en Archivo>Nuevo Proyecto y sigue las indicaciones para crear un nuevo directorio de proyecto. Luego coloca la carpeta de datos en este directorio (y asegúrate de descomprimirla).

También crea subdirectorios para ‘imágenes’ y ‘scripts’, aquí es donde guardaremos imágenes y scripts de R respectivamente.

Para proyectos complejos, querrás crear otras subcarpetas para mantenerte organizado. Por ejemplo, a menudo tengo una carpeta para documentos y una carpeta para salidas de modelos.

Ahora, cada vez que quieras trabajar en este proyecto, haz doble clic en el archivo .Rproj. Abrirá Rstudio en este directorio de trabajo (¡así que nunca necesitarás usar setwd de nuevo!).

Con el proyecto abierto, crea un nuevo script de R (Archivo>Nuevo Archivo>Script de R) y guárdalo en la carpeta de scripts. Ahora escribiremos nuestro código R en este script.

Solo verifica que estás en el directorio de trabajo de nivel superior escribiendo getwd() en la consola. Si abriste el archivo .RProj, ya estarás allí. Si no, navega allí (por ejemplo, con Sesión>Establecer Directorio de Trabajo).

Ejemplo de estructura de directorio para el proyecto

Mantenerse organizado es crítico para cualquier proyecto complejo de R. Aquí hay un ejemplo de cómo organizo mis archivos.

project/
├── data/
│   ├── raw/                # Raw data files
├── data-cleaned/            # Cleaned data files
├── scripts/
│   ├── 1_data_wrangling.R  # Script for data wrangling
│   ├── 2_dataviz.R         # Script for data visualization
│   ├── 3_modelling.R       # Script for modeling
│   └── helpers.R           # Helper functions
├── outputs/
│   ├── figures/            # Generated figures and plots
│   ├── models/             # Saved model objects
│   └── tables/             # Generated tables
├── images/                 # Images used in the project
├── docs/                   # Documentation or reports
│   └── report.Rmd          # R Markdown report
├── README.md               # Project description and instructions
└── project.Rproj           # RStudio project file

1. Recopilar los datos

¿De dónde obtengo datos para modelos?

Ya deberías tener los datos para este curso, si no descárgalos de github.

Los datos para modelos pueden provenir de muchas fuentes. La variable de respuesta es típicamente especies, presencia, abundancia o biomasa. Estos datos pueden provenir de tus propias encuestas, las de colaboradores, repositorios en línea, agencias gubernamentales o ONGs. Ejemplos de fuentes de datos incluyen encuestas pesqueras y programas de ciencia ciudadana.

También querrás algunas covariables, que usaremos para predecir las especies. Estas pueden estar incluidas con los datos, pero a menudo se obtienen de otras fuentes. Por ejemplo, la ‘Red Australiana de Datos Oceánicos’ tiene una gran cantidad de diferentes fuentes de datos oceánicos disponibles para investigadores de forma gratuita. En mi estado, el gobierno de Queensland también tiene repositorios de conjuntos de datos geográficos que se pueden descargar. Otros estados y países tienen servicios similares.

Hay paquetes especializados de R para acceder a datos ambientales, como rnaturalearth. El rLandsat puede descargar directamente datos de Landsat del USGS y ayudarte a procesarlos.

Otras fuentes de datos comunes incluyen UNEP y Zonas Económicas Exclusivas.

Si sabes qué datos quieres, intenta buscar en la web ‘r landsat’ o ‘r packages landsat’ para ver si hay un paquete disponible para facilitar la tarea de obtener los datos en un formato compatible con R.

Si no, intenta buscar el tipo de datos y ‘descargar’ o ‘shapefile’ (para datos espaciales).

Otra estrategia es mirar artículos revisados por pares que realicen análisis similares y ver qué fuentes de datos citan.

Para este curso ya te he proporcionado los datos, así que comencemos a leerlos en R.

Cargar los datos

library(tidyverse)
dat <- read_csv("data-cleaned/fish-coral-cover-sites.csv")

dat
summary(dat)

Explorando datos

Prefiero el paquete ggplot2 para graficar. Es fácil cambiar cómo se ven los gráficos con este paquete.

ggplot(dat) +
  aes(x = secchi, y = pres.topa) +
  geom_point() +
  stat_smooth()

Esto hace un gráfico de dispersión x-y. El ggplot(dat) declara el dataframe del cual extraeremos variables y también crea la página para el gráfico. El aes significa estética para el gráfico, aquí lo usamos para declarar un eje x que es la variable Abundancia de topa (pres.topa). Luego geom_points declara un objeto de geometría que decide cómo se trazan las estéticas en la página.

También hemos añadido una línea de tendencia con stat_smooth().

Hagamos también un diagrama de caja:

ggplot(dat) +
  aes(x = logged, y = pres.topa) +
  geom_boxplot()

Vemos que los topa son más abundantes en sitios no talados.

ggplot(dat) +
  aes(x = secchi, y = CB_cover) +
  geom_point() +
  stat_smooth()

dat <- dat |>
  mutate(CB_cover_percent = 100*CB_cover/n_pts)

Este código está creando una nueva columna CB_cover_percent en el dataframe dat. El operador pipe %>% se usa para pasar el dataframe dat a la función mutate, que añade una nueva columna al dataframe. La nueva columna se calcula como 100 * CB_cover / n_pts, lo que convierte la cobertura de coral a un porcentaje.

ggplot(dat) +
  aes(x = secchi, y = CB_cover_percent) +
  geom_point() +
  stat_smooth()

ggplot(dat) +
  aes(x = dist_to_logging_km, y = CB_cover_percent) +
  geom_point() +
  stat_smooth()

Explorando correlaciones en los datos

icol <- sapply(dat, is.numeric)
pairs(dat[,icol])
round(cor(dat[,icol]),2)

Que también podría resumirse como un corrplot:

icol <- sapply(dat, is.numeric)
library(corrplot)
corrplot(cor(dat[,icol]))

2. Visualización de datos

Yo haría esto en un nuevo script y lo llamaría algo como 2_dataviz.R, para que el orden de los scripts sea claro.

library(tidyverse)
dat <- read_csv("data-cleaned/fish-coral-cover-sites.csv")

Como verás, ya hemos hecho mucha visualización de datos, pero ahora haremos productos terminados para usar en un informe o publicación.

Combinando múltiples gráficos

patchwork es una forma divertida de hacer esta tarea

library(patchwork)

Primero guardamos cada ggplot como un objeto, que puede imprimirse bajo demanda:

g1 <- ggplot(dat) +
  aes(x = dist_to_logging_km, y = secchi) +
  geom_point() +
  stat_smooth() +
  xlab("Distancia a estanques de troncos (km)") +
  ylab("Profundidad Secchi (m)")  +
  theme_classic()
g1

Este código establece el tema predeterminado para todos los gráficos subsiguientes como theme_classic(), que es un tema limpio y simple.

theme_set(theme_classic())

g2 <- ggplot(dat) +
  aes(x = dist_to_logging_km, y = CB_cover) +
  geom_point() +
  stat_smooth() +
  xlab("Distancia a estanques de troncos (m)") +
  ylab("Cobertura de coral (%)")

g3 <- ggplot(dat) +
  aes(x = CB_cover, y = pres.topa) +
  geom_point() +
    stat_smooth() +
  xlab("Cobertura de coral (%)") +
  ylab("Abundancia de topa")

Ahora usa +, / y () para organizar tus gráficos en un patchwork:

gall <- g1 + g2 + g3
gall

Este código organiza los gráficos en una sola fila con los anchos especificados y agrega etiquetas a cada gráfico.

gall <- g1 + g2 + g3 +
  plot_layout(nrow = 1, widths = c(1, 0.5, 0.5)) +
  plot_annotation(tag_levels = "A")
gall

Finalmente, puedes guardar la imagen:

ggsave("outputs/plot1.png", gall, width = 8, height = 3)

Interludio: IA generativa y github copilot

Tomaremos un breve descanso del curso para ver cómo la IA generativa puede ayudarnos a hacer gráficos. Chris demostrará la herramienta github copilot.

3. Datos espaciales y creación de mapas

El paquete sf

El paquete sf se utiliza para manejar datos espaciales en R. Proporciona una interfaz simple y consistente para trabajar con datos espaciales, incluida la lectura, escritura y manipulación de objetos espaciales.

Hemos elegido la columna "pres.topa" para la escala de colores.

library(tmap)
library(sf)

#Leer tierra desde un geopackage
land <- st_read("data-cleaned/land.gpkg")

## Reading layer `land' from data source 
##   `/Users/cjbrown0/Code/R-courses/predictive-ecological-models/data-cleaned/land.gpkg' 
##   using driver `GPKG'
## Simple feature collection with 1 feature and 1 field
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 409995 ymin: 9142285 xmax: 455925 ymax: 9183985
## Projected CRS: WGS 84 / UTM zone 57S

#Leer estanques de troncos desde un geopackage
logponds <- st_read("data-cleaned/logponds.gpkg")

## Reading layer `logponds' from data source 
##   `/Users/cjbrown0/Code/R-courses/predictive-ecological-models/data-cleaned/logponds.gpkg' 
##   using driver `GPKG'
## Simple feature collection with 25 features and 2 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 425584 ymin: 9143367 xmax: 455466.5 ymax: 9163796
## Projected CRS: WGS 84 / UTM zone 57S

#Leer sdat desde un geopackage
sdat <- st_read("data-cleaned/spatial-fish-coral-cover-sites.gpkg")

## Reading layer `spatial-fish-coral-cover-sites' from data source 
##   `/Users/cjbrown0/Code/R-courses/predictive-ecological-models/data-cleaned/spatial-fish-coral-cover-sites.gpkg' 
##   using driver `GPKG'
## Simple feature collection with 49 features and 12 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 410979.7 ymin: 9143260 xmax: 454932.1 ymax: 9182978
## Projected CRS: WGS 84 / UTM zone 57S

El paquete sf proporciona una interfaz simple y consistente para manejar datos espaciales en R. Te permite leer, escribir y manipular objetos espaciales, facilitando el trabajo con datos espaciales en tus análisis.

El paquete tmap

Me gusta usar tmap para mapas.

tmap funciona muy parecido a ggplot2 en que construimos y agregamos capas. En este caso tenemos solo una capa, de sdat2. Declaramos la capa con tm_shape() (en este caso sdat2), luego el tipo de gráfico con el siguiente comando.

Aquí estamos usando tm_symbols para trazar puntos de las coordenadas. Otras opciones son tm_polygons, tm_dots y muchas otras que veremos más adelante.

Cargar nuestros datos espaciales

library(tmap)
tm_shape(sdat) +
  tm_symbols(col = "pres.topa", size = 0.2)

También podemos guardar capas de tmap, por ejemplo, la capa para la tierra será útil más adelante:

tland <- tm_shape(land) +
  tm_fill()

Ahora podemos superponer la tierra y los sitios de encuesta:

my_map <- tland +
  tm_shape(sdat) +
  tm_symbols(col = "pres.topa", size = 0.2) +
  tm_scale_bar(position = c("right", "top"))
  my_map

Lee más sobre tmap en el paquete.

Intenta hacer un mapa de una de las otras variables como secchi o CB_cover.

Para cambiar la escala de colores del mapa, puedes usar la función tm_symbols con el argumento palette. Por ejemplo:

tm_shape(sdat) +
  tm_symbols(col = "secchi", size = 0.2, palette = "Blues")

Para explorar opciones de paleta de colores:

tmaptools::palette_explorer()

tmap_save(my_map, filename = "outputs/map-pres-topa.png", width = 8, height = 6)

4. Modelado predictivo

Podrías querer comenzar otro script aquí y llamarlo 3_modelling.R.

Eligiendo un modelo

El modelado ecológico es un campo de investigación enorme y hay muchas opciones para elegir un modelo. Clasifico los modelos predictivos en dos tipos: enfoques basados en modelos y enfoques basados en aprendizaje automático (algoritmos). En este curso veremos enfoques basados en modelos, pero primero una breve descripción general.

Enfoques basados en aprendizaje automático (algoritmos)

Estos enfoques ganaron popularidad en los años 2000 con algoritmos como árboles de regresión potenciados (consulta los artículos de Jane Elith) y maxent. Son muy flexibles en términos de datos que pueden ajustar y pueden tener un excelente poder predictivo. Pero por el lado negativo, son “cajas negras” y puede ser difícil para los humanos interpretar por qué hacen buenas predicciones.

Enfoques basados en modelos

Los enfoques basados en modelos ajustan modelos a los datos. Prefiero usar enfoques basados en modelos porque me permiten interpretar los resultados (por ejemplo, ¿cómo responde una especie al aumento de temperatura?). Los enfoques de modelos también son mejores para estimar y representar límites de incertidumbre que los enfoques de aprendizaje automático. Finalmente, los enfoques basados en modelos nos permiten tratar directamente con la autocorrelación espacial. Tratar con la autocorrelación espacial es esencial para la mayoría de los conjuntos de datos ecológicos espaciales.

Otro problema que surge comúnmente (que no trataré aquí) es la detección imperfecta de especies. Los enfoques basados en modelos pueden manejar esto.

Enfoques híbridos

Siempre hay un área gris en cualquier dicotomía. Los desarrollos de modelos predictivos más recientes están utilizando redes neuronales (aprendizaje automático) para ajustar modelos complejos multiespecies y basados en procesos.

Consulta el paquete R mistnet para un ejemplo. Otro desarrollo interesante son los ‘modelos jerárquicos neuronales’ de Joseph (implementados en el lenguaje Python) que aplican aprendizaje automático para entrenar modelos de procesos de comunidades multiespecies.

Comparación de algunas opciones diferentes de modelos predictivos

Approach	Type	Pros	Cons	Example R package	When I use it
Generalized linear models	model-based	Relatively simple, good starting point, allows different response distributions (e.g. counts, binary)	Only allows linear relationships, doesn’t model spatial AC	base R, glm() function	Small or simple datasets, initial exploration of data
Generalized least squares	model-based	Relatively simple, can model spatial AC	Only allows linear relationships, only for normal distributed data	nlme, gls() function	Rarely, as I barely ever find the normal distribution is useful for ecological data
Generalized additive models	model-based	Allows for non-linear relationships, very flexible model forms can be specified, can model spatial AC and different distributions	NA	mgcv package	Basically everything these days, its such a versatile approach
Bayesian hierarchical GLMs	model-based	Allows for very flexible model forms can be specified, user can develop their own models, can model spatial AC and different distributions, can use prior information, can have GAM like splines, more completely model uncertainty	Slow computation times, requires higher technical knowledge to develop own models	brms and inla	When I want to use priori information or create my own hierarchical models (when GAMs don’t cut it)
Point process models	model-based	Useful for presence-only data	Depends on R package used for implementation, but in general this data type has lower power than presence absence or count data	spatstat, can be implmented in any package, e.g. INLA*	For presence only data
Joint models	model-based	Model multiple species and their assocations, other strenghts same as as above for Bayesian hierarchical models	As above for Bayesian models, also can be complex and very slow to fit	hmsc, boral, or write your own with stan or jags	Multispecies modelling and quantifying associations among species
Boosted regression trees	machine learning	Can handle data driven flexible non-linear relationships	No uncertainty intervals, can’t handle spatial AC (but AC can be accounted for in training design)	dismo	Only for non-spatial problems
Neural network random effects	Hybrid	Flexible relationships and multiple species, efficient for large datasets	Computationally intensive, ‘black box’	mistnet	I haven’t but am keen to try!
Neural hierarchical models	Hybrid	Flexible relationships and multiple species, efficient for large datasets, can fit process based models	Computationally intensive, ‘black box’	none yet in R	I haven’t but am keen to try!

*Nota que los modelos de procesos puntuales son matemáticamente equivalentes a maxent y pueden implementarse con muchos paquetes, por ejemplo, el enfoque SPDE del paquete INLA

Algunos ejemplos

En un estudio quería saber cómo la captura de peces se relacionaba con la abundancia de peces y variables ambientales, notablemente la temperatura. Este fue un modelo espacial y temporal de captura de peces. Desarrollé un modelo jerárquico bayesiano con el paquete brms. brms me permitió desarrollar un modelo a medida de captura en relación con la biomasa de peces en el agua y las olas de calor.

Importante, también pude modelar la incertidumbre en la variable predictora de biomasa de peces (algo que los modelos no bayesianos tienen dificultad para hacer). Los resultados sugirieron que la trucha de coral era más fácil de capturar durante una ola de calor.

Las pesquerías de palangre de atún pueden tener una alta captura incidental de especies amenazadas, como tortugas y tiburones. Entender cómo la captura incidental se relaciona con variables ambientales, pero también con especies de peces objetivo, puede ayudarnos a diseñar regulaciones más efectivas para evitar la captura incidental. Usamos modelos de distribución multiespecies (también conocidos como modelos conjuntos) para estudiar la captura de atún objetivo y captura incidental en varias pesquerías de palangre. Los modelos se implementaron con el paquete boral en R.

Los resultados revelaron grupos de captura incidental que estaban asociados con ciertas prácticas de pesca y la orientación hacia ciertas especies de atún (notablemente aleta amarilla sobre patudo).

Imagen: Ordenación de captura y captura incidental en los lances de palangre de la pesquería de atún de Palau. La ordenación se realizó a partir de un modelo conjunto ‘boral’ (Brown et al. 2021)

Finalmente, para nuestro análisis del pez loro de cabeza de joroba utilicé análisis de ruta. Esto se puede hacer con el paquete piecewiseSEM. El análisis de ruta es excelente para inferir causas de cambio (pero no tan poderoso para la predicción en mi experiencia). El análisis de ruta sugirió que los peces loro de cabeza de joroba estaban en declive tanto por la pérdida de hábitat, como por los efectos directos de la calidad del agua en los propios peces.

¿Cómo decidir qué modelo usar?

Aquí hay algunas buenas preguntas para hacer

¿Estoy más interesado en la inferencia? Si es así, evita los enfoques de aprendizaje automático
¿Estoy más interesado en hacer las predicciones más precisas, especialmente a nuevas ubicaciones? Si es así, los enfoques de aprendizaje automático suelen ser mejores.
¿Es importante representar la incertidumbre en los resultados? Si es así, evita los enfoques de aprendizaje automático y prefiere los enfoques bayesianos
¿Es la AC espacial un problema? La respuesta es casi con certeza sí, consulta la tabla anterior para ver las opciones.
¿Hay relaciones no lineales? (probablemente sí) La respuesta es casi con certeza sí, usa enfoques más flexibles como aprendizaje automático o GAMs (y enfoques bayesianos).
¿Quiero implementar fórmulas de modelo a medida? Si es así, usa enfoques bayesianos
¿Cuánto tiempo puedo esperar para los cálculos del modelo? Si es <1 minuto, entonces evita los enfoques bayesianos (o usa INLA, que es muy rápido)
¿Cuál es mi capacidad técnica? Si estás empezando, usa GLMs o GAMs. Si quieres ir por lo bayesiano, un comienzo fácil es brms o hmsc para modelos multiespecies

Preparación de covariables

Ahora que hemos cubierto algunas opciones de modelado, hagamos algo de modelado.

Vamos a comenzar con modelos más simples y avanzar. En un flujo de trabajo real, podrías hacer las cosas de manera diferente. Para un análisis complejo, es bueno comenzar simple. Mientras que para un GLM sencillo, podrías comenzar con el modelo más complejo que incluya todas las covariables y trabajar hacia atrás para simplificarlo.

Primero, crearemos variables de porcentaje de cobertura

sdat <- sdat |>
  mutate(CB_cover_percent = 100*CB_cover/n_pts,
         soft_cover_percent = 100*soft_cover/n_pts)

’

5. Modelos lineales generalizados

Si necesitas aprender rápidamente los conceptos básicos de GLMs, comienza aquí y luego lee este artículo sobre funciones de enlace

Introducción a las familias de GLM

Los Modelos Lineales Generalizados (GLMs) son una generalización flexible de la regresión lineal ordinaria que permite variables de respuesta que tienen modelos de distribución de error distintos a una distribución normal. Las diferentes familias de GLMs incluyen:

Gaussiana: Utilizada para variables de respuesta continuas con una distribución normal.
Poisson: Utilizada para datos de conteo con una distribución de Poisson.
Binomial Negativa: Utilizada para datos de conteo sobredispersos.
Binomial: Utilizada para variables de respuesta binarias.

Modelando conteos de Topa con un GLM gaussiano

Exploremos cómo se comportan diferentes familias cuando las usamos para ajustar un GLM. Modelaremos conteos de topa. Nuestra pregunta estadística será preguntar si la abundancia de Topa se ve afectada por el coral ramificado (CB_cover_percent) y/o el coral blando (soft_cover_percent). Nos gustaría saber qué hábitats son importantes para Topa.

¿Qué familia de distribución elegirías para estos datos?

m1 <- glm(pres.topa ~ CB_cover_percent + soft_cover_percent,
           data = sdat, family = "gaussian")

Este modelo intenta predecir la variable de respuesta (pres.topa) basándose en la combinación lineal de los predictores (CB_cover_percent y soft_cover_percent). Estima los coeficientes (pendientes) para cada predictor para determinar cómo los cambios en los porcentajes de cobertura de coral afectan la presencia o abundancia de Topa.

Supuestos: - La variable de respuesta (pres.topa) es continua y sigue una distribución normal. - La relación entre los predictores y la respuesta es lineal. - Homocedasticidad (varianza constante de residuos) e independencia de observaciones.

Problemas potenciales: Si pres.topa no es continua o normalmente distribuida (por ejemplo, si representa conteos o proporciones), la familia “gaussiana” puede no ser apropiada. En tales casos, otras familias como “poisson” (para conteos) o “binomial” (para resultados binarios) podrían ser más adecuadas.

Podemos verificar cómo funciona el modelo trazando algunos gráficos de diagnóstico.

par(mfrow = c(2,2))
plot(m1)

(El comando par() simplemente coloca los gráficos en un panel 2x2)

Observa cómo el gráfico QQ muestra una larga marca ascendente en valores altos. Esto significa que el modelo está subprediciendo valores altos. También nota que los residuos vs ajustados no son ‘homocedásticos’ - La dispersión de puntos de datos aumenta para valores predichos más altos.

Esto es característico de los conteos. Para sitios con baja abundancia media, la varianza es baja (porque los valores están cerca de cero). Para sitios con alta abundancia media, la varianza es alta (porque los valores están más dispersos).

Para la gaussiana también hay una inconsistencia lógica adicional - puede predecir conteos negativos. Veamos qué prediciría el modelo para la abundancia de topa a través de diferentes valores de cobertura de coral. (veremos cómo hacer estos gráficos más adelante)

Observa cómo los EE (sombreado gris) están entrando en valores negativos.

Modelando conteos de Topa con un GLM de Poisson

m2 <- glm(pres.topa ~ CB_cover_percent + soft_cover_percent,
           data = sdat, family = "poisson")

par(mfrow = c(2,2))
plot(m2)

Esto parece estar funcionando mejor, pero todavía hay algunos valores extremos.

Veamos sus predicciones:

Observa cómo las predicciones se curvan hacia arriba. Esto es porque el Poisson usa un enlace logarítmico para modelar los datos en la escala logarítmica. Así que nuestros predictores ahora son multiplicativos. Esto a menudo tiene sentido para datos de conteo.

Discutiremos la función de enlace en clase, también hay una introducción suave aquí.

Sin embargo, el Poisson todavía no está capturando toda la variabilidad en nuestros datos. Probemos la binomial negativa.

Modelando conteos de Topa con un GLM binomial negativo

library(MASS)
m3 <- glm.nb(pres.topa ~ CB_cover_percent + soft_cover_percent,
              data = sdat)

par(mfrow = c(2,2))
plot(m3)

Está funcionando mejor, aunque todavía no es perfecto. Sin embargo, nunca esperas que los datos reales sean perfectos.

La NB usa un enlace logarítmico, así que nuevamente obtenemos la curva hacia arriba. Observa cómo tiene EE mucho más amplios, que representan la variabilidad adicional en los datos y, por lo tanto, la incertidumbre en el efecto medio de CB.

Eligiendo el mejor modelo

La gaussiana es claramente inapropiada. ¿Cómo elegimos entre la Poisson y la binomial negativa?

summary(m2)

## 
## Call:
## glm(formula = pres.topa ~ CB_cover_percent + soft_cover_percent, 
##     family = "poisson", data = sdat)
## 
## Coefficients:
##                    Estimate Std. Error z value Pr(>|z|)    
## (Intercept)        0.282191   0.257565   1.096    0.273    
## CB_cover_percent   0.028180   0.004749   5.933 2.97e-09 ***
## soft_cover_percent 0.006326   0.006984   0.906    0.365    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 426.61  on 48  degrees of freedom
## Residual deviance: 372.47  on 46  degrees of freedom
## AIC: 458.37
## 
## Number of Fisher Scoring iterations: 7

summary(m3)

## 
## Call:
## glm.nb(formula = pres.topa ~ CB_cover_percent + soft_cover_percent, 
##     data = sdat, init.theta = 0.3053549986, link = log)
## 
## Coefficients:
##                     Estimate Std. Error z value Pr(>|z|)   
## (Intercept)        -0.402566   0.762327  -0.528  0.59745   
## CB_cover_percent    0.048333   0.016467   2.935  0.00333 **
## soft_cover_percent  0.008458   0.020510   0.412  0.68005   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(0.3054) family taken to be 1)
## 
##     Null deviance: 51.313  on 48  degrees of freedom
## Residual deviance: 44.168  on 46  degrees of freedom
## AIC: 207.27
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  0.3054 
##           Std. Err.:  0.0859 
## 
##  2 x log-likelihood:  -199.2710

Para el Poisson podemos comparar la desviación residual vs. los grados de libertad residuales para ver que los datos están sobredispersos para un Poisson (queremos que la proporción esté cerca de 1)

Los valores bajos de theta también indican sobredispersión. Un theta bajo significa que la NB está modelando mucha dispersión.

También podemos usar el AIC para comparar modelos.

AIC(m2, m3)

##    df      AIC
## m2  3 458.3658
## m3  4 207.2709

El modelo con el AIC más bajo es el mejor.

Más detalles sobre la teoría aquí.

Interludio: IA generativa y github copilot parte 2

Tomaremos un breve descanso del curso para ver cómo la IA generativa puede ayudarnos con la estadística ecológica. En particular, Chris mostrará cómo usarla para obtener asesoramiento estadístico sólido. Necesitas tener un poco de cuidado con lo que preguntas. Las buenas preguntas te darán buenos consejos. Las preguntas vagas o mal formuladas te darán consejos inconsistentes que a veces pueden ser incorrectos.

6. Comparando el poder predictivo de las variables de hábitat

Continuemos con la NB.

Nuestra siguiente pregunta es qué variables de hábitat importan.

Hay varias formas de hacer esto.

Comparar modelos con AIC
Usar una prueba de razón de verosimilitud
Comparar coeficientes

summary(m3)

## 
## Call:
## glm.nb(formula = pres.topa ~ CB_cover_percent + soft_cover_percent, 
##     data = sdat, init.theta = 0.3053549986, link = log)
## 
## Coefficients:
##                     Estimate Std. Error z value Pr(>|z|)   
## (Intercept)        -0.402566   0.762327  -0.528  0.59745   
## CB_cover_percent    0.048333   0.016467   2.935  0.00333 **
## soft_cover_percent  0.008458   0.020510   0.412  0.68005   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(0.3054) family taken to be 1)
## 
##     Null deviance: 51.313  on 48  degrees of freedom
## Residual deviance: 44.168  on 46  degrees of freedom
## AIC: 207.27
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  0.3054 
##           Std. Err.:  0.0859 
## 
##  2 x log-likelihood:  -199.2710

CB es significativo, pero el coral blando no lo es.

m4 <- glm.nb(pres.topa ~ CB_cover_percent,
              data = sdat)

También podríamos comparar valores de AIC:

AIC(m3, m4)

##    df      AIC
## m3  4 207.2709
## m4  3 205.4198

Así que excluir el coral blando es un modelo más parsimonioso.

anova(m3, m4, test = "Chisq")

## Likelihood ratio tests of Negative Binomial Models
## 
## Response: pres.topa
##                                   Model     theta Resid. df    2 x log-lik.   Test    df
## 1                      CB_cover_percent 0.3034218        47       -199.4198             
## 2 CB_cover_percent + soft_cover_percent 0.3053550        46       -199.2709 1 vs 2     1
##    LR stat.   Pr(Chi)
## 1                    
## 2 0.1489326 0.6995575

Esta es una prueba de razón de verosimilitud. El valor p es grande (>0.05), por lo que podemos aceptar la hipótesis nula de que los dos modelos explican una cantidad equivalente de varianza.

AIC para cobertura de coral ramificado

m5 <- glm.nb(pres.topa ~ 1, data = sdat)
AIC(m4, m5)

##    df      AIC
## m4  3 205.4198
## m5  2 209.7706

Así que solo confirmamos que es mejor tener la cobertura de coral ramificado en el modelo. Por lo tanto, m4 es nuestro “mejor” modelo.

Comparando modelos por poder predictivo

Si realmente estamos enfatizando el poder predictivo sobre la inferencia, podemos usar la validación cruzada para comparar modelos. Esta es la mejor manera de evaluar el poder predictivo. No mostraré el código ya que es un poco complejo, pero en resumen lo que hacemos es:

Dividir los datos en conjuntos de entrenamiento y prueba
Ajustar el modelo al conjunto de entrenamiento
Predecir la variable de respuesta para el conjunto de prueba
Calcular el error de predicción (por ejemplo, RMSE o MAE) para el conjunto de prueba
Repetir el proceso para cada observación en el conjunto de datos (validación cruzada leave-one-out)

El MAE se define como los valores absolutos de la diferencia entre los valores predichos y observados. Valores más bajos significan que las predicciones están más cerca de los valores observados. Echemos un vistazo a esto para nuestros dos modelos:

## Modelo 3 (con coral blando) MAE: 5.106439

## Modelo 4 (sin coral blando) MAE: 4.98766

Así que incluir coral blando aumenta el error. ¿Por qué es esto?

Una nota sobre la selección de modelos

En este caso, todos nuestros enfoques de selección de modelos nos dicen que usemos m4, el modelo sin coral blando. Pero, ¿qué pasaría si no estuvieran de acuerdo?

Discutiremos diferentes filosofías de selección de modelos en clase.

7. Haciendo predicciones a través del espacio

Gráfico predictivo con visreg

Así que hemos decidido sobre nuestro mejor modelo. Es uno con cobertura de coral ramificado y sin coral blando. Usamos la distribución NB. Veamos los efectos predichos de la cobertura de coral ramificado en la abundancia de topa.

Usaremos el paquete visreg para hacer esto.

par(mfrow = c(1,1))
library(visreg)
visreg(m4)

Esto muestra intervalos de confianza del 95% para los valores en la escala logarítmica. Discutiremos cómo interpretar esto en clase.

Aquí está el mismo modelo pero con predicciones en la escala de respuesta (es decir, número de topa por encuesta):

par(mfrow = c(1,1))
visreg(m4, "CB_cover_percent",
       gg = TRUE,
       scale = "response",
       xlab = "Cobertura de coral ramificado (%)",
       ylab = "Abundancia de topa",
       main = "Abundancia de topa predicha por cobertura de coral ramificado")

Mapa de predicciones

Una vez que tenemos nuestro mejor modelo, es sencillo hacer sus predicciones para nuestros sitios. Estas predicciones nos muestran la abundancia media esperada en un sitio dado.

sdat$predicted <- predict(m4, newdata = sdat, type = "response")

Luego podemos simplemente trazar esas predicciones en un mapa como antes:

tm_glm <- tm_shape(land) +
  tm_fill() +
  tm_shape(sdat) +
  tm_symbols(col = "predicted", size = 0.2) +
  tm_scale_bar(position = c("right", "top"))
tm_glm

Prediciendo en toda la región (interpolación)

Si quisiéramos, también podríamos usar el modelo para predecir en toda la región, de modo que podríamos hacer un mapa como este:

No cubriré el código aquí, ya que se vuelve complejo, sin embargo, aquí hay un ejemplo para este conjunto de datos en otro de mis cursos

Para hacer eso, necesitaríamos conocer nuestra covariable (cobertura CB) en todas las ubicaciones. Pero no hemos encuestado eso en todas partes.

Así que para hacer este mapa predictivo anterior, en su lugar utilicé la distancia a los estanques de troncos, que es una covariable disponible en todas partes y también está relacionada con la abundancia de topa. Luego puse la distancia a los estanques de troncos en una cuadrícula. Luego usé la misma función de predicción para predecir a esa cuadrícula. Finalmente, podemos hacer el mapa que ves arriba.

8. Modelos aditivos generalizados para relaciones no lineales

Como ejercicio final, veremos otra herramienta de modelado: ‘GAMs’

Mi favorito personal. Tan flexible, tan rápido, muy conveniente.

No puedes pasar por alto Modelos Aditivos Generalizados: Una Introducción con R, Segunda Edición si necesitas aprender sobre estos modelos.

Son realmente útiles en ecología porque permiten ajustes no lineales flexibles a los datos.

library(mgcv)
m1_gam <- gam(pres.topa ~ s(CB_cover_percent),
              family = "nb",
              data = sdat)

Se parece a nuestro GLM, sin embargo, agregamos la función s() para indicar que queremos ajustar un modelo no lineal.

Necesitamos hacer los mismos pasos que antes, así que gráficos de diagnóstico. También podríamos hacer los mismos pasos de verificación que hicimos antes.

par(mfrow = c(2,2))
gam.check(m1_gam)

## 
## Method: REML   Optimizer: outer newton
## full convergence after 3 iterations.
## Gradient range [-2.894977e-08,7.659403e-08]
## (score 98.87076 & scale 1).
## Hessian positive definite, eigenvalue range [0.5037614,10.84249].
## Model rank =  10 / 10 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##                       k'  edf k-index p-value
## s(CB_cover_percent) 9.00 2.29    0.88     0.8

Una característica agradable del AIC es que podemos comparar cualquier dos modelos siempre que tengamos los mismos datos de respuesta. Así que podemos comparar diferentes modelos, modelos con diferentes familias y modelos con diferentes covariables. Estas reglas no se aplican a las pruebas LR.

AIC(m1_gam, m4)

##              df      AIC
## m1_gam 4.822655 202.5790
## m4     3.000000 205.4198

Así que se prefiere el GAM sobre el GLM, lo que sugiere cierta no linealidad.

summary(m1_gam)

## 
## Family: Negative Binomial(0.337) 
## Link function: log 
## 
## Formula:
## pres.topa ~ s(CB_cover_percent)
## 
## Parametric coefficients:
##             Estimate Std. Error z value Pr(>|z|)   
## (Intercept)   0.7904     0.2778   2.845  0.00443 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                       edf Ref.df Chi.sq p-value   
## s(CB_cover_percent) 2.291  2.823  15.11 0.00187 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.132   Deviance explained = 25.9%
## -REML = 98.871  Scale est. = 1         n = 49

El valor edf es los grados efectivos de libertad. Esta es una medida de cuán flexible es el modelo. Valores >1 también indican no linealidad.

(gam.check(m1_gam) es una buena idea para verificar los residuos también).

Predicciones con el GAM

Podemos hacer los mismos pasos predictivos que antes con nuestro GAM para compararlo con el GLM.

g1_gam <- visreg(m1_gam, "CB_cover_percent",
       gg = TRUE,
       scale = "response",
       xlab = "Cobertura de coral ramificado (%)",
       ylab = "Abundancia de topa",
       main = "Abundancia de topa predicha por cobertura de coral ramificado")

g1_glm <- visreg(m4, "CB_cover_percent",
       gg = TRUE,
       scale = "response",
       xlab = "Cobertura de coral ramificado (%)",
       ylab = "Abundancia de topa",
       main = "Abundancia de topa predicha por cobertura de coral ramificado")

(g1_gam + g1_glm) & ylim(0, 50)

También podemos predecir a través del espacio como antes:

sdat$predicted_gam <- predict(m1_gam, newdata = sdat, type = "response")
tm_gam <- tm_shape(land) +
  tm_fill() +
  tm_shape(sdat) +
  tm_symbols(col = "predicted_gam", size = 0.2) +
  tm_scale_bar(position = c("right", "top"))

  tm_both <- tmap_arrange(tm_glm, tm_gam, ncol = 2)
  tm_both

Nota la diferencia en las escalas de color. El GAM está prediciendo abundancias más bajas en alta cobertura de coral, debido a su ajuste no lineal.

9. Conclusiones y cómo obtener ayuda

Hemos descubierto cómo la tala está impactando al pez loro de cabeza de joroba. Si quisiéramos pulir más este análisis, podríamos intersectar nuestros mapas de abundancia predicha con mapas de arrecifes (ver archivos incluidos en la carpeta de datos) y estimar cuánto hábitat de arrecife se ve afectado por la tala.

También podríamos repetir el modelado, pero usando la cobertura de coral, para hacer inferencias sobre la importancia de eso para los reclutas de topa. Ver también la sección de bonificación a continuación sobre confusión.

Obteniendo ayuda

Escribir código es 80% buscar la respuesta en Google (sin atribución)

Si vas a ser un usuario exitoso de R, entonces necesitas volverte bueno encontrando ayuda para lidiar con errores. El aforismo anterior es ampliamente suscrito por programadores profesionales. R no es una excepción. Si buscas en la web un problema, como ‘ggplot remove legend’, comúnmente obtendrás una respuesta bastante decente en Stack Overflow o un sitio similar. Probablemente usé Google decenas de veces para escribir estas notas del curso (nunca puedo recordar cómo poner los símbolos de grados en los gráficos, por ejemplo).

Si la respuesta aún no existe allí, regístrate en Stack Overflow y pregunta tú mismo (pero pasa un poco de tiempo buscando, ¡nadie quiere ser etiquetado por duplicar una pregunta existente!).

Actualizaré este aforismo:

Escribir código es 90% saber qué preguntarle a github copilot

Sin embargo, siempre ten cuidado con la IA como hemos discutido. Tiene un rendimiento asombroso (mejor que los humanos) en problemas matemáticos no resueltos. Pero las evaluaciones están encontrando que incluso los modelos más recientes son mediocres para estadísticas. Tampoco está optimizado para el programa R, por lo que todavía hay mucho con lo que no puede ayudarte (o te dará respuestas defectuosas).

Actualmente estoy trabajando en pautas para modelado ecológico con IA generativa, incluidos chatGPT y github copilot. Así que mantente atento a mi página web (https://www.seascapemodels.org/) para más información sobre eso.

Otra buena idea es encontrar un grupo de apoyo local. La codificación en R es una experiencia emocional, la frustración es común, pero la euforia de encontrar una solución puede ayudarnos a persistir. Tener otras personas que te ayuden, o incluso solo escuchen tus frustraciones, es una gran ayuda para tu motivación para seguir aprendiendo R.

PERO no te limites a búsquedas web. A menudo, la mejor manera de encontrar un método o datos es leer la literatura y ver qué han usado otras personas.

Libros de R y material web

Hay muchos buenos libros disponibles (demasiados para elegir, de hecho). Para el contenido que cubrimos hoy, algunos buenos recursos son:

R para Ciencia de Datos, principalmente para manipulación de datos
R Graphics Cookbook, para ggplot y gratis en la web
Geocomputación en R, gratis en la web
Si quieres aprender nuevos trucos, o mantenerte al día con los últimos paquetes, el agregador de blogs R-Bloggers tiene un feed continuo de blogs de R de todo el mundo y todas las disciplinas, incluido el blog de Chris

Si prefieres tener una guía impresa, otra táctica es buscar en la web tu paquete favorito y ‘cheatsheet’.

Recursos para modelado de distribución de especies

Chris a menudo se refiere a Modelos de Efectos Mixtos y Extensiones en Ecología con R. Grandes ejemplos para ayudarte a comenzar en el modelado también
El texto clásico para GAMs es Modelos Aditivos Generalizados: Una Introducción con R. Tiene un tratamiento muy técnico de los GAMs, pero también una extensa sección con ejemplos prácticos, incluidos los espaciales.
Aquí hay una introducción a algunos modelos de aprendizaje automático
Sitio para Modelado de distribución multiespecies con HMSC.

10. Material adicional

Esta es la segunda vez que imparto este curso. También hice uno anterior sobre modelado de distribución de especies que tiene más información sobre teoría y algunas otras consideraciones de modelado como autocorrelación espacial e interpolación. Así que si quieres profundizar, échale un vistazo

Validación cruzada

Aquí está el código para la validación cruzada, en caso de que estés interesado.

# Implementando validación cruzada leave-one-out (LOOCV)
n <- nrow(sdat)
predictions_m3 <- numeric(n)
predictions_m4 <- numeric(n)
observed <- sdat$pres.topa

# Bucle LOOCV
for(i in 1:n) {
  # Crear conjuntos de entrenamiento dejando fuera una observación
  train_data <- sdat[-i, ]
  test_data <- sdat[i, ]
  
  # Ajustar modelos en datos de entrenamiento
  m3_cv <- try(glm.nb(pres.topa ~ CB_cover_percent + soft_cover_percent, data = train_data), silent = TRUE)
  m4_cv <- try(glm.nb(pres.topa ~ CB_cover_percent, data = train_data), silent = TRUE)
  
  # Hacer predicciones en la observación dejada fuera
  if(!inherits(m3_cv, "try-error")) {
    predictions_m3[i] <- predict(m3_cv, newdata = test_data, type = "response")
  } else {
    predictions_m3[i] <- NA
  }
  
  if(!inherits(m4_cv, "try-error")) {
    predictions_m4[i] <- predict(m4_cv, newdata = test_data, type = "response")
  } else {
    predictions_m4[i] <- NA
  }
}

# Calcular errores de predicción
m3_rmse <- sqrt(mean((predictions_m3 - observed)^2, na.rm = TRUE))
m4_rmse <- sqrt(mean((predictions_m4 - observed)^2, na.rm = TRUE))

# Calcular error absoluto medio
m3_mae <- mean(abs(predictions_m3 - observed), na.rm = TRUE)
m4_mae <- mean(abs(predictions_m4 - observed), na.rm = TRUE)

# Imprimir resultados
cat("Modelo 3 (con coral blando) MAE:", m3_mae, "\n")
cat("Modelo 4 (sin coral blando) MAE:", m4_mae, "\n")

# Visualizar predicciones vs observados
cv_results <- data.frame(
  Observation = 1:n,
  Observed = observed,
  Model3_pred = predictions_m3,
  Model4_pred = predictions_m4
)

# Graficar observados vs predichos
ggplot(cv_results, aes(x = Observed)) +
  geom_point(aes(y = Model3_pred, color = "Modelo con coral blando"), alpha = 0.7) +
  geom_point(aes(y = Model4_pred, color = "Modelo sin coral blando"), alpha = 0.7) +
  geom_abline(intercept = 0, slope = 1, linetype = "dashed") +
  labs(
    x = "Abundancia de Topa Observada",
    y = "Abundancia de Topa Predicha",
    title = "LOOCV: Valores Observados vs Predichos",
    color = "Modelo"
  ) +
  theme_classic()

Modelado ecológico predictivo con R

CJ Brown (c.j.brown@utas.edu.au)

23 Abr 2025

Introducción

Métodos que cubriremos

Lo que necesitarás

Caso de estudio: Pez loro de cabeza de joroba, ‘Topa’ en las Islas Salomón

Planificando nuestro proyecto

1. Recopilar los datos

2. Manipulación de datos

3. Visualización de datos

4. Modelado

5. Visualización de modelos

Comencemos

Configurando tu carpeta de proyecto

Ejemplo de estructura de directorio para el proyecto

1. Recopilar los datos

¿De dónde obtengo datos para modelos?

Cargar los datos

Explorando datos

Explorando correlaciones en los datos

2. Visualización de datos

Combinando múltiples gráficos

Interludio: IA generativa y github copilot

3. Datos espaciales y creación de mapas

El paquete sf

El paquete tmap

Cargar nuestros datos espaciales

4. Modelado predictivo

Eligiendo un modelo

Enfoques basados en aprendizaje automático (algoritmos)

Enfoques basados en modelos

Enfoques híbridos

Comparación de algunas opciones diferentes de modelos predictivos

Algunos ejemplos

¿Cómo decidir qué modelo usar?

Preparación de covariables

5. Modelos lineales generalizados

Introducción a las familias de GLM

Modelando conteos de Topa con un GLM gaussiano

Modelando conteos de Topa con un GLM de Poisson

Modelando conteos de Topa con un GLM binomial negativo

Eligiendo el mejor modelo

Interludio: IA generativa y github copilot parte 2

6. Comparando el poder predictivo de las variables de hábitat

AIC para cobertura de coral ramificado

Comparando modelos por poder predictivo

Una nota sobre la selección de modelos

7. Haciendo predicciones a través del espacio

Gráfico predictivo con visreg

Mapa de predicciones

Prediciendo en toda la región (interpolación)

8. Modelos aditivos generalizados para relaciones no lineales

Predicciones con el GAM

9. Conclusiones y cómo obtener ayuda

Obteniendo ayuda

Libros de R y material web

Recursos para modelado de distribución de especies

10. Material adicional

Validación cruzada