Análisis de datos de abundancia

Práctica 3, Métodos en Biología de la Conservación

Máster en Áreas Protegidas, Recursos Naturales y Biodiversidad

Introducción

En esta práctica trabajaremos con los paquetes Distance y unmarked de R para el análisis de modelos que proporcionan estimaciones de abundancia. Analizaremos inicialmente datos procedentes de muestreos con medidas de distancia a los individuos observados, estimando densidades mediante el ajuste de las observaciones a modelos con diferentes funciones y métodos de expansión en serie. Diseñaremos los modelos, seleccionaremos los mejores e interpretaremos cuantitativamente sus resultados. Finalmente utilizaremos modelos para la estimación de densidades a partir de datos de conteos.

Preparación

Una vez iniciado R debemos cargar el archivo de datos y funciones de la asignatura (mbc.RData) accediendo al servidor con la siguiente función:

load( url( "https://webs.um.es/jfcalvo/mbc.RData" ) )

Usando ls() podemos ver los objetos cargados. Utiliza la función info para obtener información sobre ellos; por ejemplo: info( "issj" ) o info( "strix" ). Teclea info() para más información.

Necesitaremos los paquetes Distance y unmarked. Si no los tenemos instalados ejecutaremos:

install.packages( c( "Distance", "unmarked" ) )

Una vez instalados los cargaremos en memoria:

library( Distance )
library( unmarked )

Modelos de distancia con Distance

Trabajaremos inicialmente con los datos de un muestreo de distancia del argos real (Argusianus argus) en Sumatra, procedentes originalmente del sitio web acompañante del libro de Conroy y Carroll (2009). Se trata de un muestreo de 144 transectos lineales de 2200 m de longitud (considerado un ancho de banda del doble de la distancia máxima observada) en un área de estudio de 800 ha. Examinaremos en primer lugar la estructura y características de los datos:

argus

Para realizar los análisis utilizaremos la función ds del paquete Distance. Por defecto, esta función realiza ajustes para transectos lineales (argumento transect = "line"), pero puede utilizarse también para muestreos puntuales (argumento transect = "point"). Las funciones disponibles son: uniforme (key = "unif"), seminormal (key = "hn") y tasa de riesgo (key = "hr"). Los métodos de expansión en serie disponibles son: coseno (adjustment = "cos"), polinomial simple (adjustment = "pol") y polinomial de Hermite (adjustment = "her"). Ejecutaremos diversos modelos con varias de las combinaciones posibles. Por ejemplo:

ds( argus, key = "unif", adjustment = "cos" ) -> munc
ds( argus, key = "hn", adjustment = "her" ) -> mhnh
ds( argus, key = "hr", adjustment = "pol" ) -> mhrp

La selección de modelos podemos realizarla con la función summarize_ds_models:

summarize_ds_models( munc, mhnh, mhrp, output = "plain" )

Interpretaremos los resultados del mejor modelo, teniendo en cuenta la superficie total del área de estudio considerada (800 ha).

summary( mhrp )

También se pueden realizar representaciones gráficas de las funciones de detección de los modelos:

plot( mhrp )

La función ds permite incorporar covariables para modelar la detectabilidad, pero el paquete Distance no dispone de una función para realizar model averaging. No obstante, las estimaciones de model averaging podemos realizarlas en R calculando previamente los pesos de Akaike de acuerdo con la ecuación:

\[w_{i} = \frac{e^{- 0.5\ \Delta\text{AIC}_{i}}}{\sum_{}^{}e^{- 0.5\ \Delta\text{AIC}_{i}}}\] Por ejemplo, para las densidades (individuos por hectárea) estimadas por los tres modelos propuestos anteriormente:

c( 0.09858896, 0.01844645, 0.01552881 ) -> densidades
c( 0, 14.02976, 16.57292 ) -> deltas
exp( -0.5 * deltas ) / sum( exp( -0.5 * deltas ) ) -> w
sum( densidades * w )

Para obtener el intervalo de confianza hay que calcular el error estándar de la densidad media promediada (\(\bar{D}\)), utilizando los errores estándar de la densidad estimada por cada modelo (\(SE(\hat{D_i})\), con la siguiente ecuación (Conroy y Carroll 2009), :

\[SE = \sum_i{w_i \sqrt {SE(\hat{D_i})^2+(\hat{D_i}-\bar{D}})}\] \[\bar{D} \pm 1.96 \ SE\]

Ejercicio 1

Analiza los datos de un muestreo de distancias del azulillo grande (Blue Grosbeak, Passerina caerulea), realizado en Georgia, USA (79 point transects, de 100 m de radio). Estos datos proceden también del sitio web acompañante del libro de Conroy y Carroll (2009), y están disponibles en el objeto grosbeak del archivo mbc.RData. Ejecuta diversos modelos con la función ds, seleccionando las tres key functions (con las distintas series de ajuste). Recuerda usar en este caso el argumento transect = "point".

Modelos de distancia con unmarked

La función de unmarked que realiza este tipo análisis es distsamp. Esta función permite estimar la función exponencial negativa y analizar el efecto de covariables de sitio, pero requiere datos agrupados en intervalos de distancia. Analizaremos unos datos de ejemplo correspondientes a una especie de córvido, la chara de Santa Cruz (Island Scrub-jay, Aphelocoma insularis), endémica de la Isla de Santa Cruz (California). Se trata de observaciones agrupadas en tres intervalos de distancia en 307 estaciones de censo (point transects). Los datos están incluidos como ejemplo en el paquete unmarked.

head( issj )

Para su análisis deben ser transformados en un objeto de la clase unmarkedFrameDS:

issjUF <- unmarkedFrameDS( y = as.matrix( issj[ , 1:3 ] ), siteCovs = data.frame (issj[ , 6:8 ] ), dist.breaks = c( 0, 100, 200, 300), unitsIn = "m", survey = "point" )

La función distsamp no tiene opciones para diferentes métodos de expansión en serie porque utiliza un método de estimación distinto al del paquete Distance . Por lo que respecta a la key function, por defecto utiliza el argumento keyfun = "halfnorm", pero también podemos usar "hazard", "exp" o "uniform".

distsamp( ~ 1  ~ 1, data = issjUF ) -> mdu1
summary( mdu1 )

distsamp( ~ chaparral ~ 1, data = issjUF ) -> mdu2
summary( mdu2 )

distsamp( ~ chaparral ~ chaparral, data = issjUF ) -> mdu3
summary( mdu3 )

La tabla de selección de modelos la obtenemos con:

modSel( fitList( mdu1, mdu2, mdu3 ) )

Y el model averaging lo realizamos con:

head( predict( fitList( mdu1, mdu2, mdu3 ), "state" ) )

Representaremos ahora gráficamente el resultado del modelo 3, en el que tanto la densidad como la probabilidad de detección dependen de la superficie de la covariable chaparral. Para la abundancia:

plot( issj$chaparral, exp( -3.85 + 3.90 * issj$chaparral ) )

Alternativamente, podemos obtener los valores estimados por el modelo para cada unidad de muestreo (con su intervalo de confianza al 95 %) y usar la función plot.ic (disponible en el archivo mbc.RData) para representarlos:

predict( mdu3, type = "state" ) -> mdu3s
head( mdu3s )

plot.ic( x = issj$chaparral, y = mdu3s[ , c( 1, 3, 4 ) ], xlab = "Proporción de chaparral", ylab = "Abundancia" )

La representación de la probabilidad de detección requiere el uso de funciones adicionales. Para el caso de modelos estimados con la key function seminormal hay que usar la función gxhn y proporcionar el parámetro \(\sigma\) estimado, que también varía con la covariable chaparral. Así, por ejemplo, la representación de las curvas de probabilidades de detección para proporciones de chaparral iguales a 1, 0,5 y 0, respectivamente, la obtendríamos con:

plot( function( x ) gxhn( x, sigma = exp( 5.02 - 1.06 )) , 0, 100 )
plot( function( x ) gxhn( x, sigma = exp( 5.02 - 1.06 * 0.5 )) , 0, 100, add = TRUE, col = 2 )
plot( function( x ) gxhn( x, sigma = exp( 5.02 ) ), 0, 100, add = TRUE, col = 3 )

Para representar las probabilidades de detección de modelos estimados con las key functions exponencial y tasa de riesgo (hazard rate), utilizaríamos las funciones gxexp y gxhaz, respectivamente. Estas funciones se pueden utilizar también con los modelos estimados con Distance.

Ejercicio 2

Con los datos de la chara de Santa Cruz, diseña y ejecuta el mejor de los tres modelos anteriores considerando las key functions exponencial (argumento keyfun = "exp"), tasa de riesgo (argumento keyfun = "hazard") y uniforme (argumento keyfun = "uniform"). Obtén la tabla de selección del conjunto de los 4 modelos ejecutados.

Estimación de abundancia a partir de conteos

Realizaremos los análisis con datos de un muestreo de carbonero común (Parus major) consistente en 3 conteos repetidos en 263 cuadrículas de 1 km² realizados en Suiza en 2013. Estos datos proceden de un archivo que puede descargarse de la web de ejercicios del libro de Kéry y Royle (2016), en el que se recogen más especies y más años de conteos . Para nuestra práctica, los datos seleccionados (carbonero común en 2013) están disponibles en el objeto tits del archivo mbc.RData.

head( tits )

Una vez adaptados al formato correspondiente del paquete unmarked los analizaremos utilizando la función pcount.

titsUF <- unmarkedFramePCount( y = tits[ , 1:3 ], siteCovs = data.frame( elevation = tits[ , 4 ] ) )

Ejecutaremos, por ejemplo, un modelo en el que la abundancia varíe con la altitud y la detectabilidad con la visita:

pcount( ~ obsNum ~ elevation, data = titsUF ) -> mpc1
summary( mpc1 )

Los valores estimados de abundancia y probabilidad de detección se obtienen con:

predict( mpc1, type = "state" ) -> mpc1s
head( mpc1s )
predict( mpc1, type = "det" ) -> mpc1d
head( mpc1d )

La representación del submodelo de abundancia la podemos obtener con:

plot.ic( x = tits$elevation, y = mpc1s[ , c( 1, 3, 4 ) ], xlab = "Altitud (m)", ylab = "Abundancia" )

Para la representación del submodelo de detectabilidad usaremos la función plotp.ic, disponible en el archivo mbc.RData:

plotp.ic( y = mpc1d[ , c( 1, 3, 4 ) ], xlab = "Visita", ylab = "Probabilidad de detección", grupos = 1 : 3, ylim = c( 0, 1 ) )

Ejercicio 3

Con los datos del carbonero común, diseña y ejecuta dos modelos:

• uno en el que la abundancia sea dependiente de la altitud y la probabilidad de detección sea constante;

• otro en el que la abundancia muestre una respuesta unimodal frente a la altitud [~ elevation + I( elevation ^ 2 )] y la probabilidad de detección sea dependiente de la visita (obsNum).

¿Son estos modelos mejor que el mpc1? Según el mejor modelo de los tres, ¿cuál es la densidad estimada de carboneros a 1000 metros de altitud?

Ejercicio 4

Diseña y ejecuta modelos de abundancia del cárabo común (Strix aluco) en el País Vasco, considerando diferentes combinaciones para la detectabilidad y la abundancia. Los datos proceden de un estudio de Zuberogoitia et al. (2020a,b) y están disponibles en el objeto strix del archivo mbc.RData.

Evaluación

Realiza la tarea de evaluación de la práctica publicada en el Aula Virtual.

Bibliografía

• Conroy MJ, Carroll JP (2009) Quantitative conservation of vertebrates. Wiley-Blackwell, Oxford.

• Kéry M, Royle AJ (2016) Applied Hierarchical Modeling in Ecology. Volume 1. Elsevier, Amsterdam.

• Zuberogoitia I et al. (2020) Maximizing detection probability for effective large-scale nocturnal bird monitoring. Diversity and Distributions, 26: 1034-1050.

• Zuberogoitia I et al. (2020b) Maximizing detection probability for effective large-scale nocturnal bird monitoring. Dryad Digital Repository: https://doi.org/10.5061/dryad.dncjsxkwg

Descripción de los datos

Utiliza la función info. Por ejemplo: info("issj") o info("strix").