Modelo de Captura-Recaptura
Objetivos:
- Estimar el tamaño de una población usando datos de captura-marcado-recaptura, bajo dos escenarios:
- Aplicando el modelo M0.
- Modelando el sesgo asociado a la ausencia de datos de individuos que nunca fueron capturados usando la técnica de “data augmentation”.
- Determinar el área de muestreo efectiva y la densidad poblacional de la especie de interés.
Estimación del tamaño poblacional:
La estimación del tamaño de una población cerrada (i.e., que no está expuesta a migración ni mortalidad durante el período de muestreo) es un problema clásico en ecología. Para ello, hay diferentes métodos, como por ejemplo distance sampling o captura-recaptura.
En este ejercicio vamos a trabajar con datos de monito del monte (Dromiciops gliroides) en Llao Llao. Las capturas se realizaron en dos grillas de \(25\) trampas cada una. Cada individuo capturado fue marcado con un pit-tag y registramos la fecha y ubicación (trampa) de la captura.
Carguemos los datos:
trap.coor <- read.table("https://sites.google.com/site/modelosydatos/trampas.txt",
header = T)
trap.env <- read.table("https://sites.google.com/site/modelosydatos/microambientetrampa.txt",
header = T)
cap.h <- read.table("https://sites.google.com/site/modelosydatos/cap.txt", header = T)
cap.trap <- read.table("https://sites.google.com/site/modelosydatos/capturas_spat.txt",
header = T)El objeto trap.coor Son las coordenadas X e Y de las trampas. trap.envcontiene las distintas variables que describen el ambiente en donde se encuentra cada trampa. El ambiente se describe mediante variables como el grado de conectividad, el número de quintrales, la cobertura arbórea y la arbustiva y el “relleno” (que tan densa es la vegetación alrededor de la trampa).
Los datos de cap.h representan las capturas por cada individuo (en filas) para cada fecha (columnas) con \(1\) si el individuo fue capturado o con un \(0\) cuando no fue capturado. Por ejemplo, el primer individuo fue capturado solamente el 13 de Enero del 2010 y el segundo fue capturado el 11 y el 13 de Febrero.
Lo que comúnmente se obtiene a partir de estos muestreos es la historia de captura de una muestra de n individuos, el tamaño poblacional real N es desconocido, y tenemos una cantidad N-n de individuos que no fueron capturados y pertenecen a la población. Al igual que en el caso de detección imperfecta con los datos del willow tit, vamos a usar el concepto de variable oculta para poder conectar n y N. En este caso, la variable oculta va a ser la pertenencia de un individuo a la población (la llamaremos z). Si un individuo está en la población z = 1 y si no, z = 0.
El modelo más simple en la literatura de captura-recaptura se conoce como M0, y representa una conexión entre los modelos de ocupación (como el del ejemplo del willow tit) y los modelos de estimación de tamaño poblacional (ver en Royle & Dorazio 2008, pag 185). En este modelo asumimos que todos los individuos tienen la misma probabilidad de ser capturados.
Un inconveniente con el tipo de datos de captura-recaptura es que queremos estimar el tamaño de una población en base a la historia de captura de los individuos que capturamos alguna vez pero la población incluye individuos que no capturamos nunca. Para poder tener en cuenta esta posibilidad vamos a utilizar la técnica de “Data augmentation”.
Lo que hacemos con esta técnica es, literalmente, aumentar el set de datos observados con una cantidad determinada de historias de capturas de puros ceros. Estas filas completas de 0 corresponden a “pseudo-individuos” que nunca fueron capturados. Esta cantidad de pseudo-individuos, que llamaremos nz, tiene que se lo suficientemente grande para asegurarnos de que nos estamos independizando del sesgo generado por la ausencia de estos valores.
Al igual que en el ejemplo del willow tit, el desafío es conocer cuál es la fuente que genera esos ceros: son ceros generados por el muestreo o son ceros reales? El enfoque que usamos con la variable oculta es el mismo que en el ejemplo del willow tit.
cat(file = "M0_DAUG.bug", "
model{
for(i in 1: (n + nz)){
z[i] ~ dbin(psi, 1) #Variable oculta
mu[i] <- z[i] * p
y[i] ~ dbin(mu[i], J) #capturas por individuo
}
psi ~ dunif(0, 1) # probabilidad de ser parte de la población.
p ~ dunif(0, 1) # probabilidad de captura.\t
N <- sum(z[1: (n + nz)]) # Tamaño poblacional estimado.
}")Definimos las variables que vamos a usar:
n <- dim(cap.h)[1] # ind capturados
J <- dim(cap.h)[2] # total de intentos de captura
nz <- 10 # Número de ceros (individuos no capturados nunca)
# matriz de historia de capturas
Y <- rbind(as.matrix(cap.h), matrix(0, nrow = nz, ncol = J))
y <- rowSums(Y) # total de capturas por individuo Definimos la lista de datos para pasarle a JAGS, los inits y los parámetros que queremos guardar. También definimos el número de iteraciones, thinning, burnin y cadenas a correr.
m.data <- list("J", "n", "nz", "y")
inits <- function() list(psi = runif(1), p = runif(1), z = c(numeric(n) + 1,
sample(0:1, nz, replace = TRUE)))
params <- c("N", "psi", "p")
ni <- 5000
nt <- 1
nb <- 1000
nc <- 3Ahora llamamos a JAGS para que ajuste el modelo
library(jagsUI)
M0.sim <- jags(data = m.data, inits = inits, parameters.to.save = params, model.file = "M0_DAUG.bug",
n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb)##
## Processing function input.......
##
## Done.
##
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 82
## Unobserved stochastic nodes: 84
## Total graph size: 339
##
## Initializing model
##
## Adaptive phase.....
## Adaptive phase complete
##
##
## Burn-in phase, 1000 iterations x 3 chains
##
##
## Sampling from joint posterior, 4000 iterations x 3 chains
##
##
## Calculating statistics.......
##
## Done.print(M0.sim)## JAGS output for model 'M0_DAUG.bug', generated by jagsUI.
## Estimates based on 3 chains of 5000 iterations,
## adaptation = 100 iterations (sufficient),
## burn-in = 1000 iterations and thin rate = 1,
## yielding 12000 total samples from the joint posterior.
## MCMC ran for 0.04 minutes at time 2019-03-08 06:59:41.
##
## mean sd 2.5% 50% 97.5% overlap0 f Rhat n.eff
## N 80.118 1.810 76.000 81.000 82.000 FALSE 1 1.000 12000
## psi 0.966 0.029 0.892 0.974 0.999 FALSE 1 1.000 12000
## p 0.117 0.009 0.100 0.117 0.136 FALSE 1 1.001 3642
## deviance 274.033 7.369 256.567 276.609 283.446 FALSE 1 1.000 12000
##
## Successful convergence based on Rhat values (all < 1.1).
## Rhat is the potential scale reduction factor (at convergence, Rhat=1).
## For each parameter, n.eff is a crude measure of effective sample size.
##
## overlap0 checks if 0 falls in the parameter's 95% credible interval.
## f is the proportion of the posterior with the same sign as the mean;
## i.e., our confidence that the parameter is positive or negative.
##
## DIC info: (pD = var(deviance)/2)
## pD = 27.2 and DIC = 301.191
## DIC is an estimate of expected predictive error (lower is better).Luego de asegurarnos que las Rhat son menores de \(1.1\) y que los n.eff son razonables, podemos ver que el tamaño de la población muestreada es de unos \(80\) individuos. La probabilidad de captura es baja (de un \(11\)%) pero coherente para este tipo de animal. El parámetro \(\psi\) no tiene interés aquí (¿por qué?). Lo que es importante comprobar es que hayamos incluido suficientes ceros extra durante el data-augmentation. Para revisar eso vemos que la posterior de N no esté truncada a la derecha:
op <- par(cex.lab = 1.5, font.lab = 2, cex.axis = 1.3, las = 1, bty = "n")
hist(M0.sim$sims.list$N, 30, xlab = "N", main = "", freq = FALSE)
par(op)Vemos que la posterior parece estar truncada. ¿Por qué no es bueno eso? Veamos que pasa si aumentamos el valor de nz:
Ejercicio:
- Ajustar el modelo usando un
nzde \(50\) y uno de \(200\); graficar la posterior deNen cada caso y comparar.
Estimación de la densidad poblacional:
Conocer la variación en la abundancia o la densidad de una especie es relevante para la conservación y el manejo de las especies. Por ejemplo, el tamaño poblacional es el principal factor que determina si una especie se encuentra amenazada, mientras que conocer su área de acción potencial y su densidad poblacional aportan información para la toma de decisiones para su conservación.
Un problema con las estimaciones anteriores es que nos dice qué tamaño tiene la población pero no podemos relativizar al área de muestreo. Es decir, el análisis de los datos de captura-recaptura me dice que estoy muestreando una población de alrededor de \(85\) individuos, pero a qué densidad de individuos corresponde esto? Dicho de otra manera: qué área efectiva estoy muestreando con mis trampas?
El objetivo de esta parte del ejercicio es usar la ubicación de las capturas y recapturas para estimar los tamaños de áreas de acción y eventualmente la densidad poblacional. Para hacer esto, vamos a usar una modificación de la propuesta desarrollada por Gardner, Royle y Wegan. 2009. Hierarchical models for estimating density from DNA mark-recapture studies. Ecology 90:1106-1115. En este trabajo, formulan el modelo considerando que cada individuo posee un rango de acción que se representa con una normal bivariada. Entre las cosas que tenemos que estimar es dónde se encuentra el “centro de actividad” de cada individuo y qué tamaño tiene el área de acción de los individuos. La ubicación del centro de actividad afecta la probabilidad de captura.
Ahora tenemos una variable latente que es la posición del centro de actividad de cada individuo (que vamos a llamar S), y la probabilidad de detección en cada trampa depende de la distancia de la trampa al centro de actividad. El truco es asumir que el área de acción de un individuo se puede describir como una normal bivariada (porque estamos trabajando en dos dimensiones: X e Y). Entonces tenemos que estimar la ubicación de los centros de actividad y la varianza de la normal bi-variada (asumimos que es la misma para todos los individuos)
Este modelo incluye también la técnica de data augmentation que vimos en el ejemplo anterior. Ahora vamos a llamar M al total de individuos (M = n + nz). Además tenemos: X es la matriz que contiene las coordenadas en X y en Y de cada trampa. D2 es la distancia que hay entre el centro de actividad (S: Sx,Sy), y cada trampa. K es el rango de acción estimado para cada individuo en cada trampa, que decae con la distancia.
cat(file = "Mspat.bug", "
model{
for (i in 1:M) { # total de bichos, incluyendo los ceros
for (j in 1:R) { # total de trampas
D2[i,j] <- pow(Sx[i]-x[j,1], 2) + pow(Sy[i]-x[j,2], 2)
# D2 es distancia entre el centro de actividad (Sx,Sy) y la trampa j
K[i,j] <- exp(-D2[i,j]/sigma2) # gaussian kernel
gamma[i,j] <- K[i,j]/E[i] # pr de captura en la trampa j dado que
# el individuo i es capturado
}
z[i] ~ dbern(psi)
E[i] <- sum(K[i, ]) # grado de exposición del individuo i
n[i] ~ dbin(Paug[i], J) # conexión con los datos
Paug[i] <- p[i]*z[i]
p[i] <- p0 * exp(-1/E[i]) # probabilidad de captura para el individuo i
Sx[i] ~ dunif(xlower, xupper)
Sy[i] ~ dunif(ylower, yupper)
}
for(h in 1:Ncap) { # total de individuos capturados
H[h,1] ~ dcat(gamma[H[h,2], ])
# H[h,1] indica la trampa donde fue capturado el individuo H[h,2]
}
psi~ dunif(0, 1)
p0 ~ dunif(0, 1)
sigma2 ~ dunif(0, 500)
N <- sum(z[1:M])
D <- N/( (xupper-xlower) * (yupper-ylower) )
}
")Ahora le ponemos valores a algunas variables, incluyendo el número de capturas por trampa
nz <- 250 # número de ceros (individuos no capturados nunca)
Y <- rbind(as.matrix(cap.h), matrix(0, nrow = nz, ncol = J))
n <- rowSums(Y) # capturas por individuo
nind <- dim(cap.trap)[1] # total de individuos capturados
J <- dim(cap.trap)[2] # total de días de campturas
R <- dim(trap.coor)[1] # total de trampas
M <- nz + nind
# capturas por trampa
cbyt <- numeric(R)
for (i in 1:R) {
cbyt[i] <- length(which(cap.trap == i))
}Generamos una variable H donde ponemos el número de trampa en la primera columna y el ID del individuo capturado en la segunda columna. También contamos el número total de capturas Ncap.
H <- NULL #c(NA,NA)
for (i in 1:nind) {
tmp <- which(cap.trap[i, ] > 0)
for (j in 1:length(tmp)) {
H <- rbind(H, c(cap.trap[i, tmp[j]], i))
}
}
Ncap <- dim(H)[1] # toal de capturas Finalmente definimos los límites para las posibles ubicaciones de los centros de actividad y una matriz de coordenadas para las trampas
xlower <- -1.5
xupper <- 1.5
ylower <- -1.5
yupper <- 1.5
x <- as.matrix(trap.coor/100)Ahora lo de siempre; armamos la lista de datos para JAGS, la función de inits, y la lista de parámetros para guardar. Luego llamamos a JAGS (esta vez va a tardar un rato en correr).
datos <- list("xlower", "xupper", "ylower", "yupper", "x", "M", "R", "J", "Ncap",
"n", "H")
inits <- function() list(sigma2 = runif(1, 0, 50), psi = runif(1), p0 = runif(1),
z = c(numeric(nind) + 1, sample(0:1, nz, replace = TRUE)))
params <- c("sigma2", "N", "psi", "p0", "D")
Mspat.sim <- jags(data = datos, inits, parameters.to.save = params, model.file = "Mspat.bug",
n.burnin = 1000, n.chains = 3, n.iter = 2000)##
## Processing function input.......
##
## Done.
##
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 471
## Unobserved stochastic nodes: 969
## Total graph size: 182101
##
## Initializing model
##
## Adaptive phase.....
## Adaptive phase complete
##
##
## Burn-in phase, 1000 iterations x 3 chains
##
##
## Sampling from joint posterior, 1000 iterations x 3 chains
##
##
## Calculating statistics.......
##
## Done.print(Mspat.sim)## JAGS output for model 'Mspat.bug', generated by jagsUI.
## Estimates based on 3 chains of 2000 iterations,
## adaptation = 100 iterations (sufficient),
## burn-in = 1000 iterations and thin rate = 1,
## yielding 3000 total samples from the joint posterior.
## MCMC ran for 5.669 minutes at time 2019-03-08 06:59:45.
##
## mean sd 2.5% 50% 97.5% overlap0 f Rhat n.eff
## sigma2 0.148 0.019 0.118 0.146 0.192 FALSE 1 1.010 2548
## N 183.976 23.804 140.000 182.000 235.000 FALSE 1 1.004 664
## psi 0.571 0.079 0.429 0.566 0.738 FALSE 1 1.002 1084
## p0 0.174 0.021 0.138 0.173 0.218 FALSE 1 1.002 1085
## D 20.442 2.645 15.556 20.222 26.111 FALSE 1 1.004 664
## deviance 1116.023 27.437 1063.024 1115.342 1171.122 FALSE 1 1.000 2767
##
## Successful convergence based on Rhat values (all < 1.1).
## Rhat is the potential scale reduction factor (at convergence, Rhat=1).
## For each parameter, n.eff is a crude measure of effective sample size.
##
## overlap0 checks if 0 falls in the parameter's 95% credible interval.
## f is the proportion of the posterior with the same sign as the mean;
## i.e., our confidence that the parameter is positive or negative.
##
## DIC info: (pD = var(deviance)/2)
## pD = 376.4 and DIC = 1492.386
## DIC is an estimate of expected predictive error (lower is better).Si bien los Rhat son aceptables, vemos que el n.eff es medio bajo y habría que hacer más iteraciones. De todas formas vemos que la densidad poblacional es del orden de 20 individuos por hectárea.
op = par(mar = c(4, 4, 1, 1) + 0.1, las = 1, cex.lab = 1.5, cex.axis = 1.3)
hist(Mspat.sim$sims.list$D, freq = FALSE, xlab = "Individuos por ha", main = "")
par(op)En este modelo, sigma2 corresponde al cuadrado del parámetro de escala de una Weibull asumiendo forma = \(2\). Podemos estimar el tamaño de áreas de acción como:
r95 <- qweibull(0.95, shape = 2, scale = sqrt(Mspat.sim$mean$sigma2))
pi * r95^2 # home range in Ha## [1] 1.396598Si queremos calcular la posterior del área de acción hacemos:
r95s <- qweibull(0.95, shape = 2, scale = sqrt(Mspat.sim$sims.list$sigma2))
aaccs <- pi * r95s^2 # home range in Ha
op <- par(cex.lab = 1.5, font.lab = 2, cex.axis = 1.3, las = 1, bty = "n")
plot(density(aaccs), main = "", lwd = 3, xlab = "Home range (ha)", ylab = "density")
par(op)Pregunta:
- ¿Pusimos suficientes ceros extra?
sessionInfo()## R version 3.5.2 (2018-12-20)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS Mojave 10.14.3
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] jagsUI_1.5.0 lattice_0.20-38 knitr_1.20
##
## loaded via a namespace (and not attached):
## [1] Rcpp_1.0.0 bookdown_0.8 digest_0.6.18
## [4] grid_3.5.2 formatR_1.5 magrittr_1.5
## [7] coda_0.19-2 evaluate_0.12 stringi_1.2.4
## [10] rjags_4-8 rmarkdown_1.10.15 tools_3.5.2
## [13] stringr_1.3.1 parallel_3.5.2 xfun_0.4
## [16] yaml_2.2.0 compiler_3.5.2 htmltools_0.3.6Juan Manuel Morales. 6 de Septiembre del 2015. Última actualización: 2019-03-08