Este tutorial está basado principalmente en el capítulo 6 del libro Applied Hierarchical Modeling in Ecology Analysis of distribution, abundance and species richness in R and BUGS de Marc Kéry y J. Andrew Royle (Kéry and Royle 2016). Existen otros buenos tutoriales como este de M.E. Colvin o estos del Cornell Lab of Ornitology que también han sido consultados.

Introducción

La abundancia o el número de individuos de una población es un parámetro clave para conocer el estado de la misma y poder realizar una correcta gestión (Krebs 2009). Sin embargo, estimar el número total de individuos que componenen una población, \(N\), no es tarea sencilla. Para aproximarse a ese número se pueden diseñar muestreos que nos permitan obtener un índice de abundancia \(I\), que puede definirse como “cantidades que reflejan variaciones temporales o espaciales del tamaño de las poblaciones sin conocer realmente su verdadero tamaño” (Tellería 2012). De forma general, se considera que \[\begin{equation} \tag{Eq 1}\label{eq1} I = p * N \end{equation}\]

donde \(N\) es el número total de individuos de una población, \(I\) el índice de abundancia (un conteo de individuos que hayamos obtenido durante el muestreo, por ejemplo) y \(p\) sería la detectabilidad o capturabilidad, la eficacia que tenemos para detectar individuos durante los muestreos. Si logramos mantener constante esa \(p\), el índice puede ser útil para realizar comparaciones de abundancia en el espacio-tiempo. Sin embargo, muchas veces la detectabilidad está condicionada por covariables que impiden mantenerla constante (muestreos realizados en diferentes épocas del año, lugares con más o menos visibilidad, etc.). Conocer esta \(p\) y sus variaciones nos permitiría estimar el número absoluto de individuos \(N\), lo cual puede ser esencial en especies en críticamente amenazadas. Algunas aproximaciones como la captura-marcaje-recaptura se basan en la identificación individual de los individuos para estimar esta detectabilidad \(p\) y poder así obtener abundancias absolutas (Efford 2004), pero esta metodología suele ser costosa de implementar, ya que requiere una caprura y un manejo de los individuos de estudio que no es siempre posible.

Los N-mixture models, modelos de mezclas o modelo de conteos (Royle 2004) nos permiten estimar la abundancia de una especie a partir de muestreos repetidos. Estos modelos tienen una estructura jerárquica (Kéry and Royle 2016) y estiman a la vez la abundancia y la detectabilidad/capturabilidad de la especie de estudio a partir de conteos repetidos \(t\) ocasiones en \(n\) sitios diferentes.

En el siguiente tutorial se pretende comprender las generalidades de los N-mixture models a partir de la simulación de dos sencillos escenarios utilizando el lenguaje R y el paquete unmarked (Fiske and Chandler 2011).

N-mixture models: Generalidades

Pongamos que salimos al campo a buscar ciervos (Cervus elaphus).

Foto: Patricia Barroso

De forma general, se podría afirmar que si detectamos un individuo de ciervo durante el muestreo, la especie está presente en ese lugar y en ese momento (asumimos que conocemos bien la especie y que no hay posibilidad de obtener falsos positivos, esto es, confundir otras especies con individuos de ciervo, por ejemplo). Sin embargo, el hecho de no detectarla puede estar indicándonos dos cosas:

  1. Que la especie no esté presente en ese lugar porque el hábitat no es el adecuado, o la especie se haya extinguido, etc. Esto sería un verdadero negativo.
  2. Que la especie si esté presente pero no hayamos sido capaces de detectarla (hemos ido a muestrear al medio día y solo tiene actividad crepuscular, por ejemplo). Esto sería un falso negativo.

Los N-mixture models son capaces de tener en cuenta la posibilidad de cometer falsos negativos cuando muestreamos varias ocasiones los mismos sitios. Tienen en cuenta dos procesos:

  • Un proceso (o submodelo) que determina la abundancia de la especie. En inglés suele denominarse state process.
  • Un proceso (o submodelo) que determina la probabilidad de detección/observación de la especie. En inglés observational process.

Estos modelos suelen denominarse modelos jerárquicos ya que el proceso de detección depende en parte del proceso de abundancia. A continuación simularemos una serie de escenarios a modo de ejemplo para entender un poco mejor los N-mixture models.

Simulación de una población y un muestreo SIN covariables

Con ayuda de R, vamos a simular una población de ciervos que se distribuye en un área de 100 km2 (100 cuadrados de 1 km de lado)

set.seed(1) # ajustamos una "semilla"" para controlar la repetibilidad 
# Con el paquete "raster" creamos una cuadrícula de 10x10
library(raster)
sarea <- raster(nrows = 10, ncols = 10, xmn = 0, xmx = 10, ymn = 0, ymx = 10)

Este área de momento está vacía, por lo que tendremos que “llenarla” de ciervos virtuales (proceso de abundancia, o state process). Podemos utilizar una distribución de Poisson para “esparcir” nuestros ciervos simulados en el área de estudio. Esta distrubución es una de las más comunmente utilizadas en el estudio de abundancias. La distribución de Poisson maneja la frecuencia con la que un evento ocurre en un intervalo específico. Consta de un único parámetro llamado lambda \(\lambda\) que determina el número de eventos que ocurren normalmente en ese intervalo.



En nuestro caso, el intervalo una unidad espacial, una cuadrícula, mientras que el evento será la presencia de un ciervo. El parámetro \(\lambda\) será la “abundancia esperada”, esto es, la abundancia media por cuadrícula. El número y la distribución de los ciervos virtuales que colocaremos en nuestro área de estudio vendrá determinado por una distribución de Poisson con \(\lambda = 4\): \[\begin{equation} \tag{Eq 2}\label{eq2} N_{i} \sim Poisson(\lambda) \end{equation}\]

siendo \(N_{i}\) el número total de individuos en la celda \(i\).

Como podemos ver en el histograma de la izquierda, los valores más probables serán 3 o 4 ciervos, aunque variarán entre cero y 14 aproximadamente (aunque este útimo valor será muy poco probable).

# Seleccionamos el lambda deseado
lambda <- 4

# Generamos 100 números aleatorios obtenidos de una distribucion
# de Poisson con una lambda = 4
sarea[] <- rpois(100, lambda)  

# Graficamos nuestro área de estudio
plot(sarea)



Si hacemos la suma de todas las celdas de nuestro área de estudio, podremos saber el número total de individuos que hemos simulado (\(N\)):

# Calculamos el número total de ciervos virtuales creados (N)
sum(sarea[])
## [1] 368

Podemos también hacer la media de individuos en cada celda, que debería aproximarse a \(\lambda = 4\):

mean(sarea[])
## [1] 3.68


A continuación tenemos que simular nuestro muestreo (proceso de detectabilidad o observational process). Éste consistirá en 15 cuadrículas aleatorias de nuestro área de estudio a las cuales realizaremos 4 visitas (ocasiones). Durante los muestreos intentaremos detectar el máximo número posible de ciervos, pero desgraciadamente sabemos que tenemos una eficacia o probabilidad de detección del 40% para todo nuestro área de estudio. Esto significa que, en general, seremos capaces de detectar menos de la mitad de ciervos en cada muestreo. Para simular esta detectabilidad podemos utilizar una distribución binomial. La distribución binomial maneja la probabilidad de éxito para un número de sucesos con dos resultados posibles: éxito o fracaso (cara o cruz, par o impar, etc.). Pongamos que una celda tiene 10 ciervos. La detección de cada uno de esos 10 ciervos puede verse como 10 sucesos en los que puedo tener éxito (detecto al ciervo) o fracasar (el ciervo está, pero no lo detecto). De esta forma, tendríamos que \[\begin{equation} \tag{Eq 3}\label{eq3} y_{i,k} \sim binomial(N_{i}, p) \end{equation}\]

dónde \(y_{i,k}\) es el número de individuos detectados en la celda \(i\) en la visita/ocasión \(k\), \(N_{i}\) es el número de individuos en la celda \(i\) y \(p\) es la probabilidad de detectar cada ciervo, siendo en nuestro caso \(p = 0.4\). Recapitulando, muestrearemos 15 cuadrículas aleatorias durante 4 visitas diferentes. En cada visita tendremos una probabilidad del 40% de ver cada uno de los ciervos que se encuentren en la cuadrícula. El siguiente código simula este proceso:

# Seleccionamos 15 celdas aleatorias de nuestro area de estudio
site_ID <- sample(1:100, 15)
dataset <- data.frame(site_ID)

# Vamos a almacenar en número real de individuos de cada celda para poder hacer 
# comparaciones posteriormente
dataset$trueN <- extract(sarea,site_ID)

# Iniciamos nuestras 4 visitas
dataset$O1 <- NA
dataset$O2 <- NA
dataset$O3 <- NA
dataset$O4 <- NA

# Con un doble bucle, visitarremos 4 veces cada uno de nuestras 15 cuadrículas
for (j in 1:4){
  for (i in 1:length(site_ID)){
    # Nótese que los individuos detectados vienen determinados por una 
    # binomial con probabilidad 0.4, véase help(rbinom)
    dataset[i,j+2] <- rbinom(1, extract(sarea,site_ID[i]), 0.4)
  }
}

Una vez realizado el muestreo, veamos el resultado del mismo. Primero graficamos nuestro área de estudio con las celdas muestreadas:

plot(sarea)
points(xyFromCell(sarea,site_ID), pch = 4, cex = 1.5)

Ahora, vamos a explorar el resultado de nuestros muestreos en cada una de las visitas:

dataset
##    site_ID trueN O1 O2 O3 O4
## 1       19     7  3  4  0  4
## 2       50     4  1  1  2  0
## 3       95     3  1  1  0  2
## 4       77     6  3  3  4  1
## 5       87     5  2  3  2  1
## 6       75     3  0  3  2  3
## 7       91     5  2  2  3  2
## 8       33     4  1  1  2  0
## 9        7     2  1  2  0  0
## 10      38     5  2  3  2  2
## 11      32     2  0  1  0  2
## 12      96     2  1  1  0  0
## 13      54     7  1  5  2  3
## 14     100     3  1  1  0  1
## 15      90     5  2  2  3  2

La primera columna (site_ID) es un identificador de la cuadrícula. La segunda (trueN) indica el número real de individuos que se encuentran en cada cuadrícula de nuestro muestreo. Las 4 columnas siguientes (O1,O2… O4) indican el número de ciervos detectados en cada una de nuestras visitas/ocasiones. Como habíamos predicho, al tener un 40% de detectabilidad, casi siempre detectamos menos de la mitad de ciervos de los que realmente hay.

En el mundo real, el problema surge al no conocer esa detectabilidad. Al realizar los distintos muestreos no sabemos realmente cuál es el procentaje de individios que somos capaces de detectar. Una primera estima sería quedarnos con el valor máximo de cada una de nuestras visitas para cada uno de los 15 muestreos. Asumiendo que la población en cada celda es cerrada (sin migración, nacimientos o muerte de individuos) el numero máximo de individuos detectados en cada una de las 4 visitas podría ser un buen punto de partida para aproximar una abundancia total. Lo calculamos y lo comparamos con el número real del siguiente modo:

# Calculamos los individuos máximos detectados por visita
dataset$maxN <- apply(dataset[,3:6],1, max)

# Graficamos la relación con los individuos reales
plot(dataset$trueN, dataset$maxN, xlim = c(0,7.5), ylim = c(0,7.5), 
     xlab = "número real de individuos", ylab = "número detectado de individuos")
abline(0,1)

Como veníamos intuyendo, infraestimamos el número real de ciervos. De hecho, como sabemos que tenemos un 40% de eficacia de detección, podríamos adelantar que, de media, seguramente estaremos detectando aproximadamente un 40% de ciervos de los que hay en realidad. Vamos a comprobarlo:

# Calculamos la media de individuos vistos en las 4 visitas
dataset$meanN <- apply(dataset[,3:6],1, mean)

#Sumamos las medias de los 15 sitios y lo comparamos con los valores reales
sum(dataset$meanN) / sum(dataset$trueN)
## [1] 0.3928571

Efectivamente, de media estamos viendo un 39% de los individuos que realmente hay. Para poder ajustar un modelo que fuese capaz de calcular el número total de individuos, deberíamos primero ser capaces de calcular esa detectabilidad, esto es, el parámetro \(p\) de la distribución binomial.

El modelo N-mixture más sencillo con unmarked

Como ya hemos dicho, los N-mixture models nos permiten estimar a la vez la probabilidad de detección y la abundancia de nuestra especie a partir de muestreos repetidos. Para ajustar este tipo de modelos utilizaremos la función pcount del paquete unmarked. En primer lugar, debemos preparar los datos para que unmarked pueda leerlos. Utilizaremos la función unmarkedFramePCount para indicar cuales son los 15 muestreos repetidos en 4 ocasiones

library(unmarked)
dataUM <- unmarkedFramePCount(y = dataset[,3:6])

# Echamos un vistazo al objeto dataUM
head(dataUM)
## Data frame representation of unmarkedFrame object.
##    y.1 y.2 y.3 y.4
## 1    3   4   0   4
## 2    1   1   2   0
## 3    1   1   0   2
## 4    3   3   4   1
## 5    2   3   2   1
## 6    0   3   2   3
## 7    2   2   3   2
## 8    1   1   2   0
## 9    1   2   0   0
## 10   2   3   2   2

Ahora podremos ajustar nuestro modelo con la función pcount. Vamos a ajustar el modelo más sencillo posible. En este ejemplo, ni la abundancia ni la propabilidad de detección vienen determinados por ninguna covariable, simplemente responden a un proceso de Poisson con \(\lambda = 4\) y a una distribución binomial con \(p = 0.4\) respectivamente. Lo que intentaremos averiguar con este modelo será esos dos parámetros.

# Como en este modelo no hay variables predictoras, las fórmulas de detectabilidad
# y de abundancia las dejamos como ~1 respectivamente. 
# La sintaxis de pcount sería:
# pcount(fórmula_detectabilidad fórmula_abundancia, datos, número_suficientemente_alto)
# Para más información consultar help(pcount)
m1 <- pcount(~1 ~1, data=dataUM, K=50)
summary(m1)
## 
## Call:
## pcount(formula = ~1 ~ 1, data = dataUM, K = 50)
## 
## Abundance (log-scale):
##  Estimate    SE    z  P(>|z|)
##      1.54 0.373 4.12 3.75e-05
## 
## Detection (logit-scale):
##  Estimate    SE     z P(>|z|)
##    -0.596 0.561 -1.06   0.288
## 
## AIC: 185.7303 
## Number of sites: 15
## optim convergence code: 0
## optim iterations: 22 
## Bootstrap iterations: 0

La función nos devuelve las estimaciones de los parámetros \(\lambda\) y \(p\) transformados en log y logit respectivamente. Es decir, el modelo realmente ajusta el parámetro \(log(\lambda)\) y \(logit(p)\), por lo que a estas estimas hay que hacerles la transformación inversa:

# lambda estimada:
exp(coef(m1)[1]) # log
## lam(Int) 
## 4.643975
# p estimada:
#exp(coef(m1)[2])/(1+exp(coef(m1)[2])) # log ratio o log odds = logit
plogis(coef(m1)[2]) # log ratio o log odds = logit; igual que arriba
##    p(Int) 
## 0.3552852

Como vemos, la estima de \(\lambda\) es de 4.64 frente a la real, que era de 4. Por otro lado, la estima de \(p\) es de 0.35, mientras que la \(p\) real o detectabilidad era de 0.4. No es un modelo perfecto, pero se ajusta bastante a la realidad. ¿Cómo podríamos mejorar estas estimas? Incrementando el número de sitios muestreados, aumentando el número de visitas en cada sitio, etc. En cualquier caso, podemos calcular las abundancias predichas por nuestro modelo para cada sitio visitado y ver cómo se ajusta con respecto a la abundancia real. Para ello, usaremos la función bup:

# Calculamos las predicciones de abundancias
N_m1<- bup(ranef(m1)) 

# graficamos la relación entre las abundancias reales y las predichas
plot(dataset$trueN,N_m1, xlab="Abundancia simulada", ylab="Abundancia predicha")
abline(0,1)# a 1:1 line

Vemos ahora que hemos solventado el problema de la infraestima. Este sencillo ejemplo nos sirve para explicar de forma clara el funcionamiento del modelo N-mixture. Sin embargo, la abundancia normalmente suele venir determinada por alguna covariable predictora que afecta a la distribución de las especies, como disponibilidad de recursos, refugio, etc. Por eso, a continuación simularemos otro set de datos y ajustatemos un modelo un poco más complejo.

Simulación del área de estudio y del muestreo CON covariables

Supongamos ahora que nuestra población de ciervos virtuales se distribuye en un paisaje como el de la siguiente imagen



Como vemos, hay un punto de agua en la parte inferior izquierda. Vamos a crear una capa raster que indique la distancia al punto de agua, la cual usaremos posteriormente como covariable predictora para distribuir nuestros ciervos.

# Nótese que estamos estandarizando las distancias con scale(). 
# Para no obtener distancias negativas, le sumaremos 2.28668
dwat <- scale(distanceFromPoints(sarea, c(3.5,3.5))) + 2.28668

plot(dwat)

En el ejemplo anterior, habíamos distribuido los ciervos aleatoriamente siguiendo una distribución de Poisson con \(\lambda = 4\) para todo el territorio (\ref{eq2}). \[\begin{equation} \tag{Eq 2} N \sim Poisson(\lambda) \end{equation}\] Sin embargo, ahora queremos que nuestros ciervos se distribuyan conforme a la distancia al punto de agua, de tal forma que haya más ciervos cerca del agua, y su abundancia disminuya conforme nos vamos alejando de ella (estamos en una zona árida y durante la época estival el agua escasea, actuando como un atractor de ciervos). Eso quiere decir que \(\lambda\) (abundancia promedio) debe ser alta en aquellas celdas con valores bajos de distancia al agua \(dwat\), e irá disminuyendo con forme aumente \(dwat\). Este comportamiento se puede expresar como \[\begin{equation} \tag{Eq 4}\label{eq4} log(\lambda_{i}) = \beta_{0} + \beta_{1} * dwat_{i} \end{equation}\] donde \(\lambda_{i}\) es la abundancia promedio en la celda \(i\), \(\beta_{0}\) es el intercepto del modelo lineal que predice la abundancia, \(\beta_{1}\) es la “fuerza” del efecto de la distancia al agua sobre la abundancia de ciervos (la fuerza con la que el agua atrae o repele a los ciervos), y \(dwat_{i}\) es la distancia al punto de agua desde la celda \(i\). Así, podemos intuir que la \(\lambda\) de cada celda dependerá de la variable “distancia al agua”. Por lo tanto, la ecuación general que gobernará la distribución de nuestros ciervos ahora sera: \[\begin{equation} \tag{Eq 5}\label{eq5} N \sim Poisson(exp(\beta_{0} + \beta_{1} * dwat)) \end{equation}\]

Para que se cumpla nuestra premisa de que haya menos ciervos cuanto más nos alejemos del punto de agua, al coeficiente \(\beta_{1}\) deberemos darle un valor negativo (el agua atrae a los ciervos). Estos serán nuestros valores:

  • \(\beta_{0} = 3\)
  • \(\beta_{1} = -0.9\)

Calculamos la \(\lambda\) para cada celda siguiendo nuestra fórmula y los valores elegidos para los coeficientes:

beta_0 <- 3
beta_1 <- -0.9

lambda <- exp(beta_0 + beta_1*(dwat))
plot(lambda)

Una vez tenemos el \(\lambda\) de cada celda, solo nos queda “rellenar” nuestro área de estudio siguiendo, como en el ejemplo anterior, una distribución de Poisson. Sin embargo esta vez, en vez de utilizar \(\lambda = 4\) para todo el territorio, usaremos el \(\lambda\) calculado específicamente para cada una de las celdas en función de la distancia al punto de agua:

# Creamos un raster de 10x10 celdas vacías
sarea <- raster(nrows = 10, ncols = 10, xmn = 0, xmx = 10, ymn = 0, ymx = 10)  # study area

# Usamos un bucle para rellenar cada celda con el número de ciervos correspondiente según
# el lambda de cada celda
for (i in 1:ncell(sarea)){
  sarea[i] <- rpois(1, lambda[i])
}

# Graficamos la distribución de nuestros ciervos virtuales
plot(sarea)

Ya tenemos nuestros ciervos distribuidos. Recapitulando, esta población también se distribuye siguiendo una distribución de Poisson. La diferencia con respecto al ejemplo anterior es que ahora el parámentro \(\lambda\) de esa distribución no es homogéneo para todo el área de estudio, sino que varía con respecto a una covariable o variable predictora, la distancia a un punto de agua en nuestro caso. Esta situación es mucho más realista que la anterior. Podemos ver cual es el número total de individuos de nuestra población, \(N\).

sum(sarea[])
## [1] 352

Hemos simulado un total de 352 ciervos repartidos por nuestro área de estudio. Este número es el que intentaremos obtener posteriormente cuando ajustemos nuestro modelo. Para comprobar la relación entre la abundancia de ciervos y la distancia al agua, podemos graficar la correlación entre número de individuos en cada celda y la distancia de estas al punto de agua:

plot(dwat[],sarea[], xlab = "distancia al punto de agua", ylab = "nº de ciervos")

Vayamos ahora a simular nuestro muestreo. Como en el ejemplo anterior, éste consistirá en 15 cuadrículas aleatorias de nuestro área de estudio a las cuales realizaremos 4 visitas en cada una (ocasiones). Sin embargo, ahora nuestra probabilidad de detección de los animales también va a venir definido por la distancia al agua. Como vimos en la foto aérea de nuestro área de estudio, cerca del agua encontramos una mayor vegetación, lo cual dificulta la detección de nuestros ciervos, ya que se pueden ocultar mejor. Por lo tanto puestra capacidad de detección disminuirá conforme nos acerquemos al agua. Como ya habíamos visto anteriormente el proceso de observación venía definido por la fórmula (\ref{eq3}) \[\begin{equation} \tag{Eq 3} y_{i,k} \sim binomial(N_{i}, p) \end{equation}\] siendo la probabilidad de detectar un ciervo \(p\) un parámetro constante en todo nuestro área de estudio (\(p = 0.4\)) en el ejemplo anterior. Sin embargo ahora, ese parámetro \(p\) va a variar con respecto a nuestra covariable predictora de la siguiente manera: \[\begin{equation} \tag{Eq 6}\label{eq6} logit(p)= \gamma_{0} + \gamma_{1} * dwat \end{equation}\]

donde, de forma similar al modelo de \(\lambda\), ahora \(\gamma_{0}\) y \(\gamma_{1}\) son respectivamente el intercepto y la fuerza con la que la distancia al agua afecta a la probabilidad de detección de los ciervos. Pongamos que los valores son los siguientes:

  • \(\gamma_{0} = 0.2\)
  • \(\gamma_{1} = 0.8\)

Fijémnos que \(\gamma_{1}\) es positivo. Esto querrá decir que “a mayor distancia del punto de agua, mayor probabilidad de detectar un ciervo”. Vamos a calcular ahora ese parámetro \(p\) para todo nuestro área de estudio al igual que hicimos con \(\lambda\).

# logit(p) = gamma_0 + gamma_1 * (dwat)
gamma_0 <- 0.2
gamma_1 <- 0.8

# Para despejar logit() de la izquierda, pasamos a la derecha todo este cociente
p <- exp(gamma_0 + gamma_1 * dwat) / (1 + exp(gamma_0 + gamma_1 * dwat))

# Graficamos la probabilidad de detección/observación
plot(p)

Como vemos, la probabilidad de detección lejos del punto de agua llega casi al 100%. Sin embargo, cerca del punto de agua esta detectabilidad desciente hasta el 50% aproximadamente. Podemos graficar la relación entre distancia al punto de agua y la probabilidad de detección:

plot(dwat[], p[], xlab = "distancia al punto de agua", ylab = "detectabilidad", ylim = c(0.4,1))

A continuación podemos simular nuestros 15 puntos aleatorios de muestreo. Al igual que en el ejemplo anterior, visitaremos 4 veces cada punto, y aplicaremos la detectabilidad \(p\) correspondiente a cada celda.

site_ID <- sample(1:100, 15)

plot(sarea)
points(xyFromCell(sarea,site_ID), pch = 4, cex = 1.5)

dataset <- data.frame(site_ID)
dataset$trueN <- extract(sarea,site_ID)
dataset$O1 <- NA
dataset$O2 <- NA
dataset$O3 <- NA
dataset$O4 <- NA

for (j in 1:4){
  for (i in 1:length(site_ID)){
    # Nótese que utilizamos rbinom() con p = al obtenido del raster "p"
    dataset[i,j+2] <- rbinom(1, extract(sarea,site_ID[i]), extract(p,site_ID[i]))
  }
}

dataset
##    site_ID trueN O1 O2 O3 O4
## 1        5     0  0  0  0  0
## 2       79     1  1  1  1  1
## 3       66     7  6  4  5  5
## 4       87     3  3  2  3  2
## 5       27     3  3  2  3  1
## 6       95     3  3  3  1  2
## 7       88     3  3  3  2  2
## 8       91     3  3  3  2  3
## 9       37     1  1  0  1  1
## 10      71     4  3  3  4  2
## 11      65    13 10  7 10  6
## 12      72     6  5  5  4  4
## 13      24     1  1  1  1  1
## 14      89     3  2  3  2  3
## 15      32     3  3  3  2  2

Una vez obtenido nuestro conteo de ciervos, necesitamos también obtener esta vez los valores de nuestra covariable predictora para cada uno de los puntos de muestreo:

dataset$dwat <- extract(dwat, site_ID)

Modelo N-mixture con una covariable

Ya tenemos todos los elementos necesarios para ajustar un modelo N-mixture con ayuda de la función pcount de unmarked. Tenemos los conteos de ciervos en cada uno de nuestros 15 sitios durante las 4 ocasiones (visitas). Además, tenemos las medidas de una covariable que pensamos que puede ser una buena predictora de la abundancia de ciervos y, a la vez, de la detectabilidad de los mismos. Ahora debemos agrupar todos estos datos en un objeto entendible por el paquete unmarked como ya hicimos en el ejemplo anterior utilizando la función unmarkedFramePCount()

dataUM <- unmarkedFramePCount(y = dataset[,3:6], siteCovs=data.frame(dwat=dataset$dwat))

Finalmente estamos listos para ajustar un modelo en el que tengamos en cuenta nuestra covariable predictora. En este caso, tanto el modelo de abundancia como el de detectabilidad van a estar determinados por la misma covariable, la distancia al punto de agua \(dwat\), así que utilizaremos la función pcount de este modo:

# Tanto el proceso de detección como el de abundancia dependen de dwat, por tanto
# utilizaremos ~dwat ~dwat en la fórmula. Para más información sobre la 
# sintaxis de la función consultar help(pcount)
m2 <- pcount(~dwat ~dwat, data=dataUM, K=150)

# Vemos los resultados del modelo
summary(m2)
## 
## Call:
## pcount(formula = ~dwat ~ dwat, data = dataUM, K = 150)
## 
## Abundance (log-scale):
##             Estimate    SE     z  P(>|z|)
## (Intercept)    2.861 0.392  7.31 2.76e-13
## dwat          -0.872 0.207 -4.22 2.47e-05
## 
## Detection (logit-scale):
##             Estimate    SE     z P(>|z|)
## (Intercept)    0.697 0.858 0.812   0.417
## dwat           0.349 0.399 0.874   0.382
## 
## AIC: 154.0765 
## Number of sites: 15
## optim convergence code: 0
## optim iterations: 30 
## Bootstrap iterations: 0

Una diferencia importante comparado con el ejemplo previo es que anteriormente estábamos estimando los parámetros \(\lambda\) y \(p\) de los procesos de abundancia y detectabilidad (\ref{eq2} y \ref{eq3}). Como nuestro objetivo eran esos dos parámetros, debíamos “backtransformar” esas estimas, ya que el modelo los ajustaba en función de \(log(\lambda)\) y \(logit(p)\). Sin embargo, en este nuevo modelo sabemos que esos parámetros varían en función de una covariable, \(dwat\) en nuestro caso, y por tanto nuestro objetivo ya no es conocer esos dos parámetros en sí mismos, sino conocer cómo varían con respecto a la covariable predictora, esto es, buscamos los coeficientes \(\beta_{1}\) y \(\gamma_{1}\) de los modelos que gobiernan \(\lambda\) y \(p\) respectivamente (\ref{eq4} y \ref{eq6}). Por esa razón, ahora no es necesario que “backtransformemos” las estimas de esos coeficientes, sino que debemos realizar las comparaciones directamente.

En los resultados podemos ver los dos submodelos, el de abundancia y el de detectabilidad. En el de abundancia observamos una estima del intercepto \(\beta_{0} = 2.861\), mientras que la estima del efecto de la covariable predictora es \(\beta_{1} = -0.872\). Si recordamos, estos dos parámetros los habíamos simulado en \(\beta_{0} = 3\) y \(\beta_{1} = -0.9\), por lo que estamos obteniendo unas estimas bastante precisas.

Por otro lado, en el submodelo de la detectabilidad u observación, estamos obteniendo unas estimas de \(\gamma_{0} = 0.697\) y \(\gamma_{1} = 0.349\) para el intercepto y el efecto de la covariable sobre la probabilidad de detección respectivamente. Estos parámetros fueron simulados en \(\gamma_{0} = 0.2\) y \(\gamma_{1} = 0.8\), por lo que estas estimas no estarían siendo del todo acertadas En cualquier caso, sí estamos siendo capaces de detectar correctamente el efecto de la covariable \(dwat\) en ambos submodelos, negativo para el caso de la abundancia y positivo en el caso de la detectabilidad.

Una vez ajustado el modelo, podemos graficar sus predicciones teniendo en cuenta la covariable predictora \(dwat\). Utilizaremos la función predict de la siguiente forma:

# Primero debemos juntar todas las covariables en formato raster en un objeto tipo stack
# En nuestro caso solo hay una, dwat. La incluimos y le damos el nombre correspondiente
sdwat <- stack(dwat)
names(sdwat) <- "dwat"

# Ahora podemos utilizar la función predict con este stack recien creado.
# Como nuestro modelo en realidad se compone de dos submodelos, podemos realizar las 
# predicciones tanto de la abundancia como de la detectabilidad.

# type = "state" para la abundancia
pred_state <- predict(m2, newdata=sdwat, type = "state")
# type = "det" para la detectabilidad
pred_det <- predict(m2, newdata=sdwat, type = "det")

plot(pred_state)

En el gráfico superior izquierdo (Predicted) podemos ver la distribución de la abundancia predicha, mientras que en el resto se muestran el error estándar y los límites superior e inferior de las predicciones. Podemos comparar gráficamente la abundancia predicha con la abundancia simulada:

par(mfrow=c(1,2))
plot(sarea, main = "Abundancia simulada")
plot(pred_state$Predicted, main = "Abundancia predicha")

Como vemos, la distribución predicha es mucho más homogénea que la simulada, puesto que nos hemos basado en una única covariable predictora, \(dwat\). Sin embargo, el patrón general de ambas es bastante similar, por lo que podríamos considerar que nuestro modelo nos ofrece una buena predicción.

Finalmente, podemos calcular la abundancia total \(N\) predicha por nuestro modelo si hacemos el sumatorio de todas las celdas de la predicción, y compararla con la abundancia simulada:

# Recordamos el N total simulado para nuestro área de estudio
sum(sarea[])
## [1] 352
# Calculamos ahora el sumatorio de las predicciones
sum(pred_state$Predicted[])
## [1] 340.31

Vemos que la \(N\) simulada era de 352, mientras que nuestro modelo calcula una \(N\) de 340, por lo que podríamos decir que nuestras estimas están siendo bastante precisas.

Todos los parámetros estimados por nuestro modelo podrían mejorar si se aumentara el tamaño de muestreo. ¿Qué efecto tendría aumentar el número de sitios muestreados? ¿Y aumentar el número de visitas en cada sitio? ¿Afectan por igual a las estimas de los parámetros del proceso de abundancia y de observación? Como vemos, las simulaciones son una herramienta muy útil para comprender el funcionamiento de los modelos matemáticos.

Referencias

Efford, MG. 2004. “Density estimation in live-trapping studies.” Oikos 106: 598–610. doi:10.1111/j.0030-1299.2004.13043.x.

Fiske, I, and R Chandler. 2011. “unmarked: An R Package for Fitting Hierarchical Models of Wildlife Occurrence and Abundance.” Journal of Statistical Software 43 (10): 1–23. http://www.jstatsoft.org/v43/i10/.

Kéry, M, and JA Royle. 2016. “Modeling Abundance Using Binomial N-Mixture Models.” In Applied Hierarchical Modeling in Ecology Analysis of Distribution, Abundance and Species Richness in R and Bugs, 220–312. Elsevier.

Krebs, CJ. 2009. Ecology: The Experimental Analysis of Distribution and Abundance. Pearson.

Royle, JA. 2004. “N-Mixture models for estimating population size from spatially replicated counts.” Biometrics 60: 108–15. doi:10.1111/j.0006-341X.2004.00142.x.

Tellería, JL. 2012. “Prevención: Seguimiento de Especies.” In Introducción a La Conservación de Las Especies, 145–58. Tundra Ediciones.