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
Actualizado el 22/10/2023
Profesor: José Francisco Calvo Sendín
Área de Ecología, Facultad de Biología, Universidad de Murcia
jfcalvo@um.es | https:/webs.um.es/jfcalvo
1. Introducción y objetivos
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 y multi-season. Aprenderemos a diseñar e interpretar cuantitativamente los resultados de los modelos, a interpretar tablas para la selección de modelos, y a realizar estimas de probabilidad de ocupación mediante el procedimiento de model averaging.
2. Desarrollo de la sesión
(A) 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:
Con la función ls() podemos ver los datos y funciones disponibles en mbc.RData. Teclea info() para más información.
Necesitaremos además el paquete unmarked. Si no lo tenemos instalado ejecutaremos:
install.packages( "unmarked" )Una vez instalado lo cargaremos en memoria:
library( unmarked )(B) 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:
wetaUF <- unmarkedFrameOccu( weta[ , 1:5 ] )
siteCovs( wetaUF ) <- data.frame( browsed = weta[ , 6 ] )
obsCovs( wetaUF ) <- data.frame( observer = c( t( weta[ , 7:11 ] ) ) )
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 as.factor( 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
summary( m0 )Modelo 1: \(p\) dependiente de la visita, \(\psi\) constante
occu( ~ -1 + obsNum ~ 1, data = wetaUF ) -> m1
summary( m1 )Modelo 2: \(p\) dependiente de la visita, \(\psi\) dependiente del ramoneo
occu( ~ -1 + obsNum ~ browsed, data = wetaUF ) -> m2
summary( m2 )Modelo 3: \(p\) dependiente de la visita y del observador, \(\psi\) dependiente del ramoneo
occu( ~ -1 + obsNum + as.factor( observer ) ~ browsed, data = wetaUF ) -> m3
summary( m3 )Para cada sitio o unidad de muestreo, las valores estimados de \(\psi\) (argumento state) y \(p\) (argumento det) se obtienen con la función predict. Por ejemplo:
head( predict( m0, "state" ) )
head( predict( m0, "det" ) )Comprueba 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 )Comprueba también, por ejemplo, que para el 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" ) )El procedimiento de selección de modelos se realiza con:
fitList( m0, m1, m2, m3 ) -> ms
modSel( ms )
modSel( ms )@FullY el procedimiento de model averaging se realiza con:
head( predict( ms, "state") )
head( predict( ms, "det" ) )
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.0103(C) 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).
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 covertirla 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, as.factor ) )
survey <- data.frame( matrix( rep( 1:8, each = 55 ), 55, 8 * 5 ) )
survey <- data.frame( lapply( survey, as.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 1: \(\psi\) variable anualmente, \(\gamma\) constante, \(\varepsilon\) constante, \(p\) constante
colext( psiformula = ~ 1, gammaformula = ~ 1, epsilonformula = ~ 1, pformula = ~ 1, data = nsoUF) -> mm1
summary( mm1 )
predict( mm1, type = "psi" )[ 1, ]
predict( mm1, type = "col" )[ 1, ]
predict( mm1, type = "ext" )[ 1, ]
predict( mm1, type = "det" )[ 1, ]
projected( mm1 )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 coeficiesntes 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 2: \(\psi\) variable anualmente, \(\gamma\) variable anualmente, \(\varepsilon\) constante, \(p\) constante
colext(psiformula= ~ 1, gammaformula = ~ -1 + año, epsilonformula = ~ 1, pformula = ~ 1, data = nsoUF) -> mm2
summary( mm2 )
projected( mm2 )Modelo 3: \(\psi\) variable anualmente, \(\gamma\) constante, \(\varepsilon\) variable anualmente, \(p\) constante
colext(psiformula= ~ 1, gammaformula = ~ 1, epsilonformula = ~ -1 + año, pformula = ~ 1, data = nsoUF) -> mm3
summary( mm3 )
projected( mm3 )Modelo 4: \(\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) -> mm4
summary( mm4 )
projected( mm4 )Para obtener la tabla de selección de modelos y el model averaging utilizaremos:
fitList( mm1, mm2, mm3, mm4 ) -> 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, ](D) 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, as.factor ) )
survey <- data.frame( matrix( rep( 1:3, each = 100 ), 100, 3 * 5 ) )
survey <- data.frame( lapply( survey, as.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:
for (i in 1:6) print( predict( mmm0, type="phi")[[ i ]][1, ] )Los valores estimados de las probabilidades de detección son:
for (i in 1:3) print( predict( mmm0, type="det")[[ i ]][1, ] )Podemos interpretar los resultados gráficamente representando las probabilidades de ocupación y reproducción (con sus intervalos de confianza al 95%) en función del estado previo:
rbind( c( 0.3158224, 0.2353254, 0.4091226 ), c( 0.1829915, 0.07321948, 0.3883711 ), c( 0.9510666, 0.8137157, 0.9885688) ) -> psi
psi
barplot( psi[ , 1 ], ylim = c( 0, 1 ) ) -> bp
segments( bp, psi[ , 2 ], bp, psi[ , 3 ])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)3. Evaluación
Realiza la tarea de evaluación de la práctica publicada en el Aula Virtual.
4. 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.
Descripción de los datos
Utiliza la función info. Por ejemplo: info("nso").