Datos

Para la elaboración de los modelos de distribución a partir del muestreo de campo realizado, se utilizó la información de 277 cuadrículas UTM 1×1 km. Durante los muestreos, se intentó completar cada cuadrícula con al menos un itinerario de unos 15 min. o 500 m. En total se consiguieron 286 itinerarios/transectos con una longitud total de aproximadamente 145 km.

La información recogida en cada transecto sobre la abundancia de especies y hábitats se utilizó para realizar un análisis sobre la selección de hábitat de cada especie (véase más abajo). Por otro lado, para realizar los modelos de distribución (o de probabilidad de presencia) se utilizó la información de presencia o ausencia de las especies a escala de celda UTM de 1×1 km tanto del trabajo de campo directo como de otras fuentes que pudieran determinar la presencia de especies a esa escala. Esto permitió aumentar el número de observaciones de cada especie mediante la anexión de observaciones de presencia descargadas de la plataforma eBird (eBird Basic Dataset, 2024). En concreto, se descargaron datos para los años 2020 a 2024 en los periodos entre abril-junio para toda la Sierra de Guadarrama y procesados con el paquete de R auk (Strimas-Mackey et al., 2023).

Para cada celda UTM de 1×1 km se extrajeron todas las observaciones de cada especie y se asignaron a las categorías detectada o no detectada. Además, se realizaron modelos sólo para aquellas especies que al menos se hubiesen detectado en al menos 5 celdas UTM en los transectos propios del trabajo de campo de este altas para asegurar un número mínimo de datos para la aplicación de los modelos.

 

Variables ambientales

Para los modelos se utilizaron diferentes variables que describen distintos aspectos del territorio como la topografía del terreno, el uso del suelo o la climatología. Estas variables fueron extraídas de diferentes fuentes de datos y convertidas a escala de 1km:

• Elevación media proveniente del Multi-Error-Removed Improved-Terrain Digital Elevation Model (MERIT-DEM, Yamazaki et al., (2017).
• Índice de vegetación mejorado (EVI en inglés) derivado de imágenes del satélite Landsat 8 cortesía del servicio geológico de EE.UU. y ajustando lo más posible el periodo de los muestreos con el periodo de las imágenes, es decir, primavera de 2024.
• Ocupación del suelo y porcentaje de hábitat principal clasificado por los observadores durante los muestreos y obtenido para las celdas UTM de 1×1 del Sistema de Información sobre Ocupación del Suelo de España de Alta Resolución (SIOSE-AR) para el año 2017.
• Variables bioclimáticas derivadas de valores mensuales de temperatura y precipitaciones.

 

Modelos de distribución

En este trabajo se aplicaron modelos de distribución basados en datos de presencia y ausencia para identificar las condiciones ambientales asociadas a la probabilidad de presencia de la especie en cada tipo de hábitat. Cabe destacar que estos modelos no informan directamente sobre la distribución espacial del número de individuos, aunque puede asumirse intuitivamente que una alta probabilidad de presencia podría estar asociada a una alta abundancia, esto no necesariamente es cierto. La probabilidad de presencia indica simplemente que las condiciones permiten la existencia de la especie, mientras que la abundancia depende de procesos adicionales como la disponibilidad de recursos, la competencia entre individuos de la misma especie y entre otras especies, o la calidad del hábitat.

Para realizar unos buenos modelos de abundancia que informen fehacientemente de dónde se encuentran más o menos aves, se necesita un muestreo repetido de la zona a estudiar a lo largo de los años, cosa que en el presente trabajo no se pudo realizar.

En ecología, es habitual que tanto las especies como las variables ambientales muestren patrones de autocorrelación espacial, lo que significa que los sitios geográficamente próximos tienden a ser más similares entre sí que aquellos más alejados (Dormann et al., 2007). Si esta autocorrelación espacial no se controla en los modelos estadísticos, puede llevar a estimaciones sesgadas. Por ello, para controlar la autocorrelación espacial debidamente, se utilizaron modelos espacialmente explícitos que incorporan efectos aleatorios espaciales usando el paquete de R sdmTMB (Anderson et al., 2024).

De acuerdo con los requisitos del modelo, se creó una malla utilizando el algoritmo k-means de agrupamiento para definir la distancia mínima permitida entre puntos en las unidades de X e Y, en km (véase Anderson et al. (2024) para más detalles).

La creación de efectos aleatorios espaciales permite incorporar en el modelo efectos espaciales que no son capturados por las variables ambientales pero que determinan la distribución de las especies (por ejemplo, si una especie es gregaria o si decide asentarse en una zona concreta por razones desconocidas).

Para los modelos de distribución o de probabilidad de presencia se utilizó una distribución binomial convirtiendo las abundancias a Presencia/Ausencia (o de 1 y 0, respectivamente). Para incorporar posibles relaciones no lineales de las variables ambientales con las distribuciones de las especies, se utilizaron Modelos Generales Aditivos (GAM en inglés).

Además, para reducir la colinealidad entre variables, se seleccionaron solo aquellas con un coeficiente de correlación de Pearson menor a 0,7, y reteniendo aquellas con mayor sentido biológico general. Todas las variables se escalaron (es decir, se centró la media) antes de incluirlas para facilitar la interpretación de los coeficientes y mejorar la convergencia del modelo (enfoque de modelo completo, (Whittingham et al., 2006)).

Debido a que cada especie tiene unos requerimientos ecológicos diferentes, se aplicaron una media de 100 modelos para cada una de ellas y se seleccionó el modelo con un BIC menor (Bayesian Information Criterion) y que, por tanto, se ajustaba mejor a los datos.
Después de ajustar el modelo, se utilizó la función sanity() de sdmTMB para verificar la validez y estabilidad de las estimaciones. Esta función realiza una serie de comprobaciones esenciales que se detallan en el paquete (Anderson et al., 2024).

Para evaluar la capacidad predictiva del mejor modelo, se ajustó con el 80 % de los datos para el entrenamiento del modelo y el 20 % para probarlo (Elith y Leathwick, 2009) y se calcularon los parámetros AUC, Sensibilidad (o Tasa de Verdaderos Positivos, TPR) y Especificidad (o Tasa de Verdaderos Negativos, TNR).

Por último, para hacer una representación de la probabilidad de presencia más detallada, se desagregaron los mapas probabilísticos de la escala 1x1km a una escala de 100×100 m.

Selección de hábitat

Aunque la información de abundancias recogida en sólo un año de campo es insuficiente para realizar estimaciones de la abundancia de especies, sí que permite realizar análisis exploratorios de la relación entre la densidad de las diferentes especies y el hábitat que ocupan. Para ello, de cada transecto se extrajo la abundancia total observada de cada especie (es decir, la suma total de individuos), y se calculó la densidad de individuos por 500 metros. De esta manera, las abundancias en transectos con diferentes longitudes son comparables.

Para estudiar la selección de hábitat de las especies, se tradujeron las diferentes categorías de hábitat recogidas por los voluntarios en el campo a categorías más generales que permiten visionar las relaciones mejor (figura 1). Posteriormente se utilizaron diagramas de cajas para observar estas relaciones.

 

Figura 1. Mapa de distribución de cada tipo de hábitat considerado en el Parque Nacional de la Sierra de Guadarrama correspondiente a la Comunidad de Madrid.