load( url( "https://webs.um.es/jfcalvo/mbc.RData" ) )
Análisis de datos de ocupación
Práctica 2, 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 el paquete unmarked de R para el análisis de modelos de ocupación. Utilizaremos diversos datos de ejemplo para realizar análisis single-season, de comunidades, multi-season y multi-state. Aprenderemos a diseñar e interpretar cuantitativamente los resultados de los modelos, a interpretar tablas para la selección de modelos, y a realizar estimaciones de probabilidad de ocupación mediante el procedimiento de model averaging.
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:
Usando ls() podemos ver los objetos cargados. Utiliza la función info para obtener información sobre ellos; por ejemplo: info( "weta" ) o info( "nso" ). Teclea info() para más información.
Necesitaremos además los paquetes unmarked y AICcmodavg. Si no los tenemos instalados ejecutaremos:
install.packages( c( "unmarked", "AICcmodavg" ) )Una vez instalados los cargaremos en memoria:
library( unmarked )
library( AICcmodavg )Modelos single-season. Ejemplo del weta
Trabajaremos inicialmente con los datos del weta gigante (Mahoenui giant weta, Deinacrida mahoenui) un ortóptero de la familia Anostostomatidae endémico de Nueva Zelanda. La presencia de esta especie fue muestreada en 72 parcelas circulares de 3 m de radio, visitadas en 5 días distintos de marzo de 2004. Como covariable de sitio se considera la existencia de ramoneo por cabras del arbusto Ulex europaeus. El ramoneo favorece el follaje de la planta, que aporta alimento y protección ante depredadores. Como covariable de muestreo se considera el efecto de 3 observadores distintos.
wetaPara trabajar con las funciones de unmarked necesitamos convertir la tabla de datos en un objeto de la clase unmarkedFrameOccu, especificando los datos de ocupación, las covariables de los sitios o unidades de muestreo y las covariables correspondientes a cada día de observación:
datos <- weta [ , 1 : 5 ]
browsed <- data.frame( browsed = weta[ , 6 ] )
observer <- data.frame( weta[ , 7 : 11 ] )
observer <- data.frame( lapply( observer, factor ) )
wetaUF <- unmarkedFrameOccu( y = datos, siteCovs = browsed, obsCovs = list( observer = observer ) )
wetaUF
summary( wetaUF )Diseñaremos y ejecutaremos diversos modelos con la función occu, con diferentes especificaciones para estimar la probabilidad de detección y la probabilidad de ocupación (en este orden). Los efectos del ramoneo se especifican con la variable browsed, los efectos del día de visita con obsNum y los del observador con observer. Con el argumento -1 el modelo no estima el intercept, por lo que obtenemos un coeficiente propio para cada día de observación. Una vez realizados los análisis procedamos a: (1) interpretar los resultados de cada modelo, (2) analizar la tabla de selección de modelos y sus parámetros, y (3) realizar el model averaging de las estimaciones.
Modelo nulo: \(p\) constante, \(\psi\) constante
occu( ~ 1 ~ 1, data = wetaUF ) -> m0
m0Modelo 1: \(p\) dependiente de la visita, \(\psi\) constante
occu( ~ -1 + obsNum ~ 1, data = wetaUF ) -> m1
m1Modelo 2: \(p\) dependiente de la visita, \(\psi\) dependiente del ramoneo
occu( ~ -1 + obsNum ~ browsed, data = wetaUF ) -> m2
m2Modelo 3: \(p\) dependiente de la visita y del observador, \(\psi\) dependiente del ramoneo
occu( ~ -1 + obsNum + observer ~ browsed, data = wetaUF ) -> m3
m3Los valores estimados de \(\psi\) (argumento state) y \(p\) (argumento det) se obtienen con la función predict. Por ejemplo:
predict( m0, "state" )
predict( m0, "det" )Observa que en el caso de \(\psi\) obtenemos 72 valores (para cada uno de los sitios o unidades de muestreo), mientras que en el caso de \(p\) obtenemos 360 (uno para cada sitio y visita: \(72 \times 5\)). Podemos comprobar que ambas probabilidades se calculan a partir de los coeficientes aplicando la transformación logística inversa mediante la función plogis.
summary( m0 )
plogis( 0.475124 )
plogis( -0.621831 )En el caso del modelo 2, las probabilidades de ocupación de sitios ramoneados y no ramoneados son, respectivamente:
summary( m2 )
plogis( -0.027318 + 1.235149 )
plogis( -0.027318 )
head( predict( m2, "state" ) )Observa también que los valores estimados de detectabilidad varían para cada visita:
head( predict( m2, "det" ) )Podemos representar los resultados usando la función plotp.ic (disponible en el archivo mbc.RData) con los valores estimados por predict. Por ejemplo, las probabilidades de ocupación:
predict( m2, "state" ) -> m2ps
plotp.ic( m2ps[, c( 1, 3 : 4 ) ], ylim = c( 0, 1 ), ylab = "Probabilidad de ocupación", grupos = c( "Ramoneado", "No ramoneado" ) )Y las probabilidades de detección:
predict( m2, "det" ) -> m2pd
plotp.ic( m2pd[, c( 1, 3 : 4 ) ], ylim = c( 0, 1 ), ylab = "Probabilidad de detección", grupos = c( "Visita 1", "Visita 2", "Visita 3", "Visita 4", "Visita 5" ) )El procedimiento de selección de modelos se realiza con:
fitList( m0, m1, m2, m3 ) -> ms
modSel( ms )Los resultados completos (y con mayor número de decimales) se obtienen con:
modSel( ms )@FullFinalmente, el procedimiento de model averaging se realiza con:
predict( ms, "state" )
predict( ms, "det" )Observa que el caso de las probabilidades de detección tenemos múltiples valores que corresponden a diferentes observadores y visitas. Podemos visualizarlos con mayor claridad si escribimos lo siguiente:
t( matrix( ( predict( ms , "det" )[ , 1]), 5, 72) )Comprueba que la probabilidad de ocupación estimada en sitios ramoneados, de acuerdo con el model averaging es:
0.7674 * 0.6180 + 0.7699 * 0.2466 + 0.6298 * 0.1251 + 0.6166 * 0.0103Pruebas de bondad de ajuste
La función mb.gof.test del paquete AICcmodavg permite realizar pruebas de bondad de ajuste de los modelos de ocupación single-season y multi-season (siguiente sección), con la estimación del factor de inflación de la varianza \(\hat{c}\). Por ejemplo:
mb.gof.test( m3, nsim = 100 )También dispone de una función que proporciona una tabla de selección de modelos basada en el AIC corregido (AICc) y en el quasi AIC (QAICc):
aictab( c( m0, m1, m2, m3 ) )
aictab( c( m0, m1, m2, m3 ), c.hat = 1.3 )Observa las diferencias entre estas tablas y la proporcionada por la función modSel de unmarked.
Modelos de ocupación de comunidades
Cuando se trabaja con múltiples especies en estudios de comunidades, adquiere especial interés la estimación de los valores de riqueza (número de especies), que debido a la detección imperfecta son siempre superiores a los valores observados. La función occuComm del paquete unmarked permite realizar este tipo de análisis. Como ejemplo utilizaremos los datos de presencia/ausencia de 60 especies de aves en itinerarios de censo de 1 km de longitud (3 visitas), realizados en 34 localizades (zonas agroforestales) de la Región de Murcia. Los datos, procedentes de una tesis doctoral realizada en la Universidad de Murcia (Zamora 2021), están disponibles en el objeto avesUF, de la clase unmarkedFrameOccuComm.
avesUFPor defecto solo se muestran los datos de la primera especie, junto con la covariable ambiental altitud y la covariable de observación (visita). También se incluye la riqueza observada de cada itinerario (riq.obs) con objeto de realizar la comparación con las estimaciones de los modelos.
Si queremos visualizar los datos de todas las especies debemos utilizar:
avesUF@ylistLa función occuComm realiza análisis de ocupación para el conjunto de las especies. Por ejemplo:
occuComm( ~ visita ~ altitud, data = avesUF ) -> mc1
mc1Para cada especie, las probabilidades de ocupación estimadas por el modelo, se obtienen con predict:
lapply( predict( mc1, "state"), head )También las probabilidades de detección:
lapply( predict( mc1, "det"), head )Si queremos obtener las estimaciones de una especie concreta podemos usar el nombre o su índice La estimación de la riqueza para cada localidad se obtiene con:
richness( mc1 )De esta forma, podemos comparar riquezas observadas y estimadas mediante el siguiente gráfico:
plot( avesUF@siteCovs$riq.obs, richness( mc1 ) )
abline( 0, 1 )Modelos multi-season. Ejemplo del búho moteado
Utilizaremos en este caso los datos de ejemplo del búho moteado (Strix occidentalis caurina) correspondendientes a un estudio realizado en California entre 1997 y 2001, en el que 55 territorios fueron visitados hasta en 8 ocasiones durante la época de reproducción para determinar la presencia de los búhos, es decir, la ocupación de cada territorio (objeto nso).
.jpg#/media/File:Northern_Spotted_Owls_(51983452420).jpg){kind=link}
nso
dim( nso )Los dos nuevos parámetros que aparecen en este tipo de modelos son \(\gamma\) (probabilidad de colonización) y \(\varepsilon\) (probabilidad de extinción local). Aun sin covariables, el mayor número de parámetros a estimar determina que el número total de modelos que podrían considerarse como candidatos sea elevado. Por tanto, en general, la investigación debe centrarse en un conjunto más reducido de modelos atendiendo a hipótesis concretas. Por ejemplo, podemos limitar el número de hipótesis sobre las probabilidades de detección asumiendo que estas son constantes (no varían entre visitas ni entre años). En consecuencia, todos los modelos que vamos a ejecutar inicialmente cuentan con una estructura fija para las probabilidades de detección (\(p\) constante).
Al igual que en la sección anterior, debemos transformar la tabla de datos original (nso) para convertirla en este caso en un objeto de la clase unmarkedMultFrame:
year <- data.frame( matrix( rep( 1997 : 2001, each = 55 ), 55, 5 ) )
year <- data.frame( lapply( year, factor ) )
survey <- data.frame( matrix( rep( 1 : 8, each = 55 ), 55, 8 * 5 ) )
survey <- data.frame( lapply( survey, factor ) )
nsoUF <- unmarkedMultFrame( y = nso, obsCovs = list( visita = survey ), yearlySiteCovs = list( año = year ), numPrimary = 5 )
nsoUF
summary( nsoUF )Procederemos a diseñar y ejecutar los siguientes modelos:
Modelo 0 (nulo): \(\psi\) variable anualmente, \(\gamma\) constante, \(\varepsilon\) constante, \(p\) constante
colext( psiformula = ~ 1, gammaformula = ~ 1, epsilonformula = ~ 1, pformula = ~ 1, data = nsoUF) -> mm0
summary( mm0 )
predict( mm0, type = "psi" )[ 1, ]
predict( mm0, type = "col" )[ 1, ]
predict( mm0, type = "ext" )[ 1, ]
predict( mm0, type = "det" )[ 1, ]
projected( mm0 )Este modelo estima la probabilidad de ocupación del primer año (1997). Los parámetros derivados (calculados apartir de los parámetros estimados) son las probabilidades de ocupación a partir del segundo año (\(\psi_{1998 - 2001}\)) y las tasas de cambio anual (\(\lambda_{t}\)):
\[\psi_{t + 1} = \psi_{t} + {\ \gamma}_{t}\left( 1 - \psi_{t} \right) - \varepsilon_{t}{\ \psi}_{t}\]
\[\lambda_{t} = \frac{\psi_{t + 1}}{\psi_{t}}\] Comprueba que \({\ \psi}_{1997}\), \(\gamma\), \(\varepsilon\) y \(p\) se calculan a partir de los coeficientes del modelo:
plogis( 0.537 )
plogis( -1.49 )
plogis( -1.73 )
plogis( -0.021 )Comprueba también que, por ejemplo, \({\ \psi}_{1998}\) y \(\lambda_{1997}\) se calculan como:
0.63116 + 0.18418 * ( 1 - 0.63116 ) - 0.15074 * 0.63116
0.6039519 / 0.63116Modelo 1: \(\psi\) variable anualmente, \(\gamma\) variable anualmente, \(\varepsilon\) constante, \(p\) constante
colext( psiformula= ~ 1, gammaformula = ~ -1 + año, epsilonformula = ~ 1, pformula = ~ 1, data = nsoUF ) -> mm1
summary( mm1 )
projected( mm1 )Comprueba que, en este caso, las tasas de colonización son diferentes cada año:
predict( mm1, type = "col" )Modelo 2: \(\psi\) variable anualmente, \(\gamma\) constante, \(\varepsilon\) variable anualmente, \(p\) constante
colext(psiformula= ~ 1, gammaformula = ~ 1, epsilonformula = ~ -1 + año, pformula = ~ 1, data = nsoUF ) -> mm2
summary( mm2 )
projected( mm2 )Modelo 3: \(\psi\) variable anualmente, \(\gamma\) variable anualmente, \(\varepsilon\) variable anualmente, \(p\) constante
colext(psiformula= ~ 1, gammaformula = ~ -1 + año, epsilonformula = ~ -1 + año, pformula = ~ 1, data = nsoUF ) -> mm3
summary( mm3 )
projected( mm3 )Para obtener la tabla de selección de modelos y el model averaging utilizaremos:
fitList( mm0, mm1, mm2, mm3 ) -> mms
modSel( mms )
predict( mms, type = "psi" )[ 1, ]
predict( mms, type = "col" )[ 1 : 5, ]
predict( mms, type = "ext" )[ 1 : 5, ]
predict( mms, type = "det" )[ 1 : 40, ]Modelos multi-season multi-state
Analizaremos unos datos simulados sobre la ocupación y reproducción de una especie en 100 territorios, visitados en 3 ocasiones durante 5 años (objeto msms). Utilizaremos la función occuMS del paquete unmarked y ejecutaremos un único modelo estimando las probabilidades condicionales de ocupación y reproducción en función del estado previo de territorio. Los estados son: no ocupado (0), ocupado sin reproducción (1), ocupado con reproducción (2). Así, por ejemplo, \(\psi^{\lbrack 0\rbrack}\) (phi[0] en R) denota la probabilidad de ocupación de un territorio en el año t + 1 dado que no estaba ocupado en el año \(t\), y \(R^{\lbrack 2\rbrack}\) (R[2] en R) denota la probabilidad de reproducción de un territorio en el año \(t + 1\) dado que estaba ocupado con reproducción en el año \(t\).
De esta forma, es posible establecer una matriz de transición (proceso Markoviano), mediante la que se definen las diferentes probabilidades de cambio de estado de los territorios entre el año \(t\) (columnas) y el año \(t + 1\) (filas):
\[\text{A} = \begin{bmatrix} 1 - \psi^{\lbrack 0\rbrack} & 1 - \psi^{\lbrack 1\rbrack} & 1 - \psi^{\lbrack 2\rbrack} \\ \psi^{\lbrack 0\rbrack}\left( 1 - R^{\lbrack 0\rbrack} \right) & \psi^{\lbrack 1\rbrack}\left( 1 - R^{\lbrack 1\rbrack} \right) & \psi^{\lbrack 2\rbrack}\left( 1 - R^{\lbrack 2\rbrack} \right) \\ \psi^{\lbrack 0\rbrack}R^{\lbrack 0\rbrack} & \psi^{\lbrack 1\rbrack}R^{\lbrack 1\rbrack} & \psi^{\lbrack 2\rbrack}R^{\lbrack 2\rbrack} \end{bmatrix}\]
Por otra parte, existen tres parámetros de detección distintos: \(p_{t, j}^{\lbrack{1}\rbrack}\), la probabilidad de detectar la ocupación en la visita \(j\) del año \(t\) en un territorio ocupado sin reproducción en el año \(t\); \(p_{t, j}^{\lbrack{2}\rbrack}\), la probabilidad de detectar la ocupación en la visita \(j\) del año \(t\) en un territorio ocupado con reproducción en el año \(t\); y \(\delta_{t,j}\), la probabilidad de detectar la reproducción en la visita \(j\) del año \(t\) en un territorio ocupado con reproducción.
En primer lugar transformaremos la tabla de datos en un objeto de clase unmarkedFrameOccuMS, añadiendo las covariables de año y visita:
msms
year <- data.frame( matrix( rep( 1 : 5, each = 100 ), 100, 5 ) )
year <- data.frame( lapply( year, factor ) )
survey <- data.frame( matrix( rep( 1 : 3, each = 100 ), 100, 3 * 5 ) )
survey <- data.frame( lapply( survey, factor ) )
msmsUF <- unmarkedFrameOccuMS( y = msms, obsCovs = list( visita = survey ), yearlySiteCovs = list( año = year ), numPrimary = 5 )
summary( msmsUF )Diseñaremos y ejecutaremos un único modelo en el que tanto la ocupación como la reproducción dependan del estado previo del territorio. Asumiremos también probabilidades de detección constantes, independientes de la visita y del año.
psiformulas <- c( "~ 1", "~ 1" ) # psi inicial, R inicial
phiformulas <- c( "~ 1", "~ 1", "~ 1", "~ 1", "~ 1", "~ 1" ) # psi[0], psi[1], psi[2], R[0], R[1], R[2]
detformulas <- c( "~ 1", "~ 1", "~ 1" ) # p[1], p[2], delta
occuMS( detformulas, psiformulas, phiformulas, data = msmsUF, parameterization = "condbinom" ) -> mmm0
summary( mmm0 )Los valores estimados de las probabilidades condicionadas de ocupación y reproducción son:
lapply( predict( mmm0, type = "phi" ), head )Los valores estimados de las probabilidades de detección son:
lapply( predict( mmm0, type = "det" ), head )Podemos interpretar los resultados gráficamente representando las probabilidades de ocupación, reproducción y detección (con sus intervalos de confianza al 95 %) en función del estado previo. Usaremos la función plotp.ic. Por ejemplo:
rbind( c( 0.3158224, 0.2353254, 0.4091226 ), c( 0.1829915, 0.07321948, 0.3883711 ), c( 0.9510666, 0.8137157, 0.9885688) ) -> psi
psi
plotp.ic( psi, ylim = c( 0, 1 ), xlab = "Estado previo", ylab = "Probabilidad de ocupación", grupos = c( "0", "1", "2" ) )Comprobaremos también que la matriz de transiciones se elabora a partir de los parámetros estimados según:
A <- matrix(NA, 3, 3)
A[ 1, ] <- 1 - c( 0.3158224, 0.1829915, 0.9510666 )
A[ 2, ] <- c( 0.3158224, 0.1829915, 0.9510666 ) * ( 1 - c( 0.1662257, 0.6848002, 0.8767159 ) )
A[ 3, ] <- c( 0.3158224, 0.1829915, 0.9510666 ) * c( 0.1662257, 0.6848002, 0.8767159 )
round ( A, 3 )
apply( A, 2, sum )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.
Fiske I, Chandler R (2011) unmarked: An R Package for Fitting Hierarchical Models of Wildlife Occurrence and Abundance. Journal of Statistical Software 43: 1-23.
Kéry M, Royle AJ (2016) Applied Hierarchical Modeling in Ecology. Volume 1. Elsevier, Amsterdam.
MacKenzie DI, Nichols JD, Royle JA, Pollock KH, Bailey LL, Hines JE (2018) Occupancy Estimation and Modeling - Inferring Patterns and Dynamics of Species Occurrence. 2nd ed. Academic Press.
Zamora JM (2021) Ecological contributions of small waterbodies to animal biodiversity in a Mediterranean semiarid region. Tesis doctoral, Universidad de Murcia.
Descripción de los datos
Utiliza la función info. Por ejemplo: info( "weta" ) o info( "nso" ).