Modelo de Ocupación Dinámico

Diego J. Lizcano

Este documento contiene un ejemplo del modelo de ocupación de una sola especie y varias temporadas. Este modelo fue descrito por MacKenzie et. al (2003) en:

MacKenzie, D. I., J. D. Nichols, J. E. Hines, M. G. Knutson, and A. B. Franklin. 2003. Estimating site occupancy, colonization, and local extinction when a species is detected imperfectly. Ecology 84:2200–2207.

image

Parámetros Adicionales

Este modelo incluye dos parámetros adicionales que representan la colonización y extinción de cada sitio. Estos dos parámetros adicionales pueden ser modelados con covariables que varían año a año.

Antes de entrar en materia con el modelo dinámico vale la pena que revisemos el modelo básico de ocupación.

Recordemos el Modelo estático

Recordemos que el modelo básico de ocupación de MacKenzie et. al (2002), también es conocido como el modelo estático de ocupación. Este modelo se aplica a una sola especie, y por lo general en una sola temporada.

image

Donde
ψ es la ocupación y p la probabilidad de detección. Con
β como el coeficiente de la regresión para las co-variables de la ocupación y
α el coeficiente de regresión para las co-variables de la detección.

Si desean conocer en detalle el modelo estático y saber más del poder y potencial de las simulaciones en ecología, les recomiendo seguir el Tutorial Ubicado en este enlace.

Recordemos la forma de la distribución Bernoulli!

La distribución de Bernoulli (o distribución dicotómica), nombrada así por el matemático suizo Jacob Bernoulli (1655-1705).

image

Esta es una distribución de probabilidad discreta, que toma valor 1 para la probabilidad de éxito (p) y valor 0 para la probabilidad de fracaso (q=1-p).

La distribución de Bernoulli es un caso particular de la distribución binomial, pero con solo dos posibles resultados (éxito o fracaso) unos y ceros.

El proceso Bernoulli es el más simple proceso aleatorio que existe! Imagínemelo como algo tan sencillo como una secuencia de lanzamientos de una moneda. Donde un solo lanzamiento es un “trial” y muchos lanzamientos componen el proceso.

Usemos el siguiente código ejecutándolo varias veces y cambiando los parámetros para entender su efecto en el resultado.

library(ggplot2)
ni<-10 # numero de datos
pi<- 0.5 # probabilidad (~proporcion de unos)
# Generemos datos con esa informacion 
daber<-data.frame(estimado=rbinom(ni, 1, pi)) 
# Grafiquemos 
library(ggplot2)
ggplot(daber, aes(x=estimado)) + 
    geom_histogram(aes(y=..density..), # Histograma y densidad 
                   binwidth=.1, # Ancho del bin
                   colour="black", fill="white") + 
        geom_vline(aes(xintercept=mean(estimado, na.rm=T)), 
          color="blue", linetype="dashed", size=1) # media en azul

Si desean conocer más detalles de la distribución Bernoulli les recomiendo la muy buena explicación de jbstatistics

Y si quieren aún más detalles visiten la clase del profesor Tsitsiklis del MIT.

Los Parámetros de los modelos

Modelo estatico:

\(z_{i} \sim Bernoulli (\psi_{i})\) Proceso Ecológico

\(y_{ij} \sim Bernoulli (z_{i} * p_{ij})\) Proceso de Observación

con sitio i durante el muestreo j

De forma Lineal

logit(Ψi) = α0 + α1xi1 + . . . + αUxiU.

logit(pij) = β0 + β1xi1 + . . . + βUxiU + βU+1yij1 + . . . + βU+V yijV.

Modelo Dinámico

\(z_{it} \sim Bernoulli (\psi_{i,t})\) Proceso Ecológico

\(y_{ijt} \sim Bernoulli (z_{i} * p_{i,tj})\) Proceso de Observación

con sitio i durante el muestreo j en el tiempo t (años)

Modelo Dinámico considerando colonización y extinción

image

\(z_{it} \sim Bernoulli (z_{i,t-1\phi it } + (1-z_{i,t-1}) \gamma_{i,t})\) Proceso Ecológico

\(y_{ijt} \sim Bernoulli (z_{it} * p_{ijt})\) Proceso de Observación

De forma Lineal

logit(Ψi,1) = α0 + α1xi1 + . . . + αUxiU.

Ψi,t|zi,t−1 = zi,t−1 × (1 −\(\epsilon_{i}\),t−1) + (1 − zi,t−1) × \(\gamma_{i}\),t−1, for t > 1.

La función Colext del paquete unmarked

image

Veamos una simulacion con la función Colext de unmarked

Modelo de Ocupación Diámico con datos simulados

Primero generamos un conjunto de datos simple y simulado con valores específicos de año específicos para los parámetros, así como especificaciones de diseño, es decir, número de sitios, años y encuestas por año.

Luego, veremos cómo ajustar un modelo de ocupación dinámico con dependencia del año en los parámetros de probabilidad de colonización, extinción y detección.

Simulando, formateando y resumiendo datos

Para simular los datos, ejecutamos el siguiente código R. Los valores reales para estos parámetros para cada año se extraen aleatoriamente de una distribución uniforme con los límites especificados.

M <- 250 # Number of sites  
J <- 3 # num secondary sample periods 
T <- 10 # num primary sample periods  
psi <- rep(NA, T) # Occupancy probability  
muZ <- z <- array(dim = c(M, T)) # Expected and realized occurrence  
y <- array(NA, dim = c(M, J, T)) # Detection histories  

set.seed(13973)  
psi[1] <- 0.4 # Initial occupancy probability  
p <- c(0.3,0.4,0.5,0.5,0.1,0.3,0.5,0.5,0.6,0.2)  
phi <- runif(n=T-1, min=0.6, max=0.8) # Survival probability (1-epsilon)
gamma <- runif(n=T-1, min=0.1, max=0.2) # Colonization probability  

# Generate latent states of occurrence  
# First year  
z[,1] <- rbinom(M, 1, psi[1]) # Initial occupancy state  
# Later years  
for(i in 1:M){ # Loop over sites 
  for(k in 2:T){ # Loop over years 
    muZ[k] <- z[i, k-1]*phi[k-1] + (1-z[i, k-1])*gamma[k-1] 
    z[i,k] <- rbinom(1, 1, muZ[k]) 
    } 
  }  # Generate detection/non-detection data  
for(i in 1:M){ 
  for(k in 1:T){ 
    prob <- z[i,k] * p[k] 
    for(j in 1:J){ 
      y[i,j,k] <- rbinom(1, 1, prob) 
      } 
    } 
  }  # Compute annual population occupancy  
for (k in 2:T){ 
  psi[k] <- psi[k-1]*phi[k-1] + (1-psi[k-1])*gamma[k-1] 
  }

Hemos generado una sola realización del sistema estocástico así definido. La Figura 1 ilustra la cuestión fundamental de la detección imperfecta: la proporción real de sitios ocupados difiere mucho de la proporción observada de sitios ocupados, y debido a que p varía entre años, los datos observados no pueden usarse como un índice válido del parámetro de interés ψi.

plot(1:T, colMeans(z), type = "b", xlab = "Year", ylab = "Proportion of sites occupied", col = "black", xlim=c(0.5, 10.5), 
      xaxp=c(1,10,9), ylim = c(0,0.6), 
      lwd = 2, lty = 1, frame.plot = FALSE, 
      las = 1, pch=16)  

psi.app <- colMeans(apply(y, c(1,3), max))  
lines(1:T, psi.app, type = "b", col = "blue", 
      lty=3, lwd = 2)  
legend(1, 0.6, c("truth", "observed"), 
      col=c("black", "blue"), lty=c(1,3), pch=c(16,1))

Analizando los datos

Para analizar este conjunto de datos con un modelo de ocupación dinámico sin marcar, primero cargamos el paquete, y luego a continuación, formateamos los datos de detección / no detección de una matriz tridimensional (como se genera) en una matriz bidimensional con M filas. Es decir, colocamos las tablas de datos anuales (los segmentos de la matriz 3-D anterior) de lado para producir un diseño “amplio” de los datos. Posteriormente, creamos una matriz que indica el año en que se muestreo cada sitio.

library(unmarked)
yy <- matrix(y, M, J*T)
year <- matrix(c('01','02','03','04','05','06','07','08','09','10'), 
               nrow(yy), T, byrow=TRUE)

Para organizar los datos en el formato requerido por colext, utilizamos la función unmarkedMultFrame. Los únicos argumentos requeridos son y, los datos de detección / no detección, y numPrimary, el número de estaciones. Los tres tipos de covariables descritos anteriormente también se pueden suministrar utilizando los argumentos siteCovs, annualSiteCovs y obsCovs. En este caso, solo usamos el segundo tipo, que debe tener M filas y T columnas.

simUMF <- unmarkedMultFrame(y = yy, 
                            yearlySiteCovs = list(year = year), 
                            numPrimary=T)

summary(simUMF)

unmarkedFrame Object

250 sites Maximum number of observations per site: 30 Mean number of observations per site: 30 Number of primary survey periods: 10 Number of secondary survey periods: 3 Sites with at least one detection: 195

Tabulation of y observations: 0 1 6430 1070

Yearly-site-level covariates: year
01 : 250
02 : 250
03 : 250
04 : 250
05 : 250
06 : 250
(Other):1000

Construcción y ajuste de modelos

Estamos listos para adaptar algunos modelos de ocupación dinámica. Ajustaremos un modelo con valores constantes para todos los parámetros (modelo nulo) y otro con dependencia total del tiempo para la probabilidad de colonización, extinción y detección.

 # Model with all constant parameters  
m0 <- colext(psiformula= ~1, 
             gammaformula = ~ 1, 
             epsilonformula = ~ 1, 
             pformula = ~ 1, 
             data = simUMF, 
             method="BFGS")

summary(m0)

Call: colext(psiformula = ~1, gammaformula = ~1, epsilonformula = ~1, pformula = ~1, data = simUMF, method = “BFGS”)

Initial (logit-scale): Estimate SE z P(>|z|) -0.813 0.158 -5.16 2.46e-07

Colonization (logit-scale): Estimate SE z P(>|z|) -1.77 0.0807 -22 2.75e-107

Extinction (logit-scale): Estimate SE z P(>|z|) -0.59 0.102 -5.79 7.04e-09

Detection (logit-scale): Estimate SE z P(>|z|) -0.0837 0.0562 -1.49 0.137

AIC: 4972.597 Number of sites: 250 optim convergence code: 0 optim iterations: 27 Bootstrap iterations: 0

El tiempo de cálculo fue de solo unos segundos. Tenga en cuenta que todos los parámetros se estimaron en la escala logit. Para volver a transformar a la escala original, simplemente podemos usar la función de logit inverso, llamada plogis en R.

Alternativamente, podemos usar backTransform, que calcula los errores estándar usando el método delta. Los intervalos de confianza también se obtienen fácilmente utilizando la función con fi n. Primero nos recordamos los nombres de los parámetros, que pueden usarse como argumentos para estas funciones.

plogis(-0.813)

[1] 0.3072516

names(m0)

[1] “psi” “col” “ext” “det”

backTransform(m0, type="psi")

Backtransformed linear combination(s) of Initial estimate(s)

Estimate SE LinComb (Intercept) 0.307 0.0335 -0.813 1

Transformation: logistic

confint(backTransform(m0, type="psi"))  # intervalos
 0.025     0.975

0.2457313 0.3765804 Luego ajustamos el modelo de ocupación dinámica con dependencia total del año en los parámetros que describen la dinámica de ocupación y también en la detección. Este es el mismo modelo bajo el cual generamos el conjunto de datos, por lo que esperaríamos estimaciones precisas.

Por defecto en R, un factor como el año en este análisis, se parametriza en términos de una intersección y efectos que representan diferencias. Esto significaría que el parámetro para el primer año es la intersección y los efectos denotarían las diferencias entre los valores de los parámetros en todos los demás años, en relación con el valor del parámetro en el primer año, que sirve como nivel de referencia. Este tratamiento o la parametrización de los efectos es útil para evaluar las diferencias. Para una presentación simple, una parametrización de medias es más práctica. Se puede especificar agregando un -1 a la fórmula para los parámetros dependientes del tiempo

 m1 <- colext(psiformula = ~1, # First-year occupancy 
              gammaformula = ~ year-1, # Colonization 
              epsilonformula = ~ year-1, # Extinction 
              pformula = ~ year-1, # Detection 
              data = simUMF)

m1

Call: colext(psiformula = ~1, gammaformula = ~year - 1, epsilonformula = ~year - 1, pformula = ~year - 1, data = simUMF)

Initial: Estimate SE z P(>|z|) -0.273 0.302 -0.906 0.365

Colonization: Estimate SE z P(>|z|) year01 -2.08 0.951 -2.19 2.86e-02 year02 -2.18 0.365 -5.96 2.52e-09 year03 -1.98 0.274 -7.23 4.88e-13 year04 -2.32 0.678 -3.42 6.37e-04 year05 -1.89 0.478 -3.95 7.78e-05 year06 -1.76 0.294 -5.97 2.44e-09 year07 -1.55 0.230 -6.73 1.75e-11 year08 -1.43 0.228 -6.29 3.19e-10 year09 -2.35 0.470 -5.00 5.64e-07

Extinction: Estimate SE z P(>|z|) year01 -1.4209 0.418 -3.401 6.72e-04 year02 -0.4808 0.239 -2.009 4.45e-02 year03 -1.2606 0.366 -3.440 5.83e-04 year04 -0.0907 0.650 -0.139 8.89e-01 year05 -0.6456 0.599 -1.078 2.81e-01 year06 -0.9586 0.378 -2.539 1.11e-02 year07 -1.2279 0.365 -3.362 7.74e-04 year08 -1.1894 0.292 -4.076 4.58e-05 year09 -0.6292 0.635 -0.991 3.22e-01

Detection: Estimate SE z P(>|z|) year01 -1.0824 0.244 -4.434 9.26e-06 year02 -0.2232 0.148 -1.508 1.32e-01 year03 0.2951 0.154 1.918 5.52e-02 year04 0.0662 0.161 0.412 6.81e-01 year05 -2.0396 0.433 -4.706 2.52e-06 year06 -0.6982 0.232 -3.005 2.66e-03 year07 0.2413 0.165 1.466 1.43e-01 year08 0.0847 0.155 0.548 5.84e-01 year09 0.6052 0.140 4.338 1.44e-05 year10 -1.1699 0.306 -3.828 1.29e-04

AIC: 4779.172

Selección de Modelos

models <- fitList(
    'psi(.)gam(.)eps(.)p(.)' = m0,
    'psi(.)gam(year)eps(year)p(year)' = m1
      )

ms <- modSel(models)
ms
                            nPars     AIC  delta AICwt cumltvWt

psi(.)gam(year)eps(year)p(year) 29 4779.17 0.00 1e+00 1.00 psi(.)gam(.)eps(.)p(.) 4 4972.60 193.43 1e-42 1.00

Predicción y Graficas

Nuevamente, todas las estimaciones se muestran en la escala logit. Las estimaciones de transformación inversa cuando hay covariables, como el año, implican un paso adicional. Específicamente, necesitamos decir sin marcar los valores de nuestra covariable en los que queremos una estimación. Esto se puede hacer usando backTransform en combinación con linearComb, aunque puede ser más fácil de usar predict.

predic le permite al usuario proporcionar un marco de datos en el que cada fila representa una combinación de valores covariables de interés. A continuación, creamos los data.frames llamados nd y cada fila representa un año. Luego solicitamos estimaciones anuales de la probabilidad de extinción, colonización y detección, y las comparamos con la “verdad”, es decir, los valores con los que simulamos el conjunto de datos. Tenga en cuenta que hay parámetros de extinción y colonización T-1 en este caso, por lo que no es necesario incluir el año “10” en nd.

Predict es mas versatil y devuelve las predicciones junto con errores estándar e intervalos de confianza. Estos se pueden usar para crear graficas. La función with se usa para simplificar el proceso de solicitud de las columnas de data.frame devueltas por predic.

# Crear nuevo data frame
nd <- data.frame(year=c('01','02','03','04','05','06','07','08','09'))
# predecir
E.ext <- predict(m1, type='ext', newdata=nd)
E.col <- predict(m1, type='col', newdata=nd) 

nd <- data.frame(year=c('01','02','03','04','05','06','07','08','09','10')) 
E.det <- predict(m1, type='det', newdata=nd)

## Graficas 
### Extinction
# op <- par(mfrow=c(3,1), mai=c(0.6, 0.6, 0.1, 0.1)) 
with(E.ext,{ # Plot for extinction probability 
  plot(1:9, Predicted, pch=1, xaxt='n', xlab='Year', 
       ylab=expression(paste('Extinction probability ( ', epsilon, ' )')), 
       ylim=c(0,1), col=4)
  
  axis(1, at=1:9, labels=nd$year[1:9]) 
  arrows(1:9, lower, 1:9, upper, code=3, angle=90, length=0.03, col=4)
  points((1:9)-0.1, 1-phi, col=1, lwd = 1, pch=16) 
  legend(7, 1, c('Parametro verdadero', 'Estimado'), col=c(1,4), pch=c(16, 1), cex=0.8) 
  }) 

### colonization
 with(E.col, { # Plot for colonization probability 
   plot(1:9, Predicted, pch=1, xaxt='n', xlab='Year', ylab=expression(paste('Colonization probability ( ', gamma, ' )')), ylim=c(0,1), col=4) 
   axis(1, at=1:9, labels=nd$year[1:9]) 
   arrows(1:9, lower, 1:9, upper, code=3, angle=90, length=0.03, col=4)
   points((1:9)-0.1, gamma, col=1, lwd = 1, pch=16) 
   legend(7, 1, c('Parameter', 'Estimate'), col=c(1,4), pch=c(16, 1), cex=0.8) 
   })

### Detection
with(E.det, { # Plot for detection probability: note 10 years 
  plot(1:10, Predicted, pch=1, xaxt='n', xlab='Year', 
       ylab=expression(paste('Detection probability ( ', p, ' )')), 
       ylim=c(0,1), col=4) 
  
  axis(1, at=1:10, labels=nd$year)
  arrows(1:10, lower, 1:10, upper, code=3, angle=90, length=0.03, col=4)
  points((1:10)-0.1, p, col=1, lwd = 1, pch=16) 
  legend(7.5, 1, c('Parameter','Estimate'), col=c(1,4), pch=c(16, 1), cex=0.8) 
  })

# par(op)

Probando la bondad del ajuste del modelo

Además de estimar la varianza de una estimación, la rutina de bootstrap paramétrica se puede usar para evaluar la bondad del ajuste. Para este propósito, una estadística de ajuste, es decir, una que compara los valores observados y esperados se evalúan utilizando el modelo ajustado original y muchos otros modelos ajustados a conjuntos de datos simulados. La simulación produce una aproximación de la distribución de la estadística de ajuste, y un valor P puede calcularse como la proporción de valores simulados mayor que el valor observado. Hosmer y col. (1997) encontraron que un estadístico χ2 funcionó razonablemente bien al evaluar la falta de ajuste para los modelos de regresión logística. No conocemos estudios que evalúen formalmente el desempeño de varias estadísticas de ajuste para modelos de ocupación dinámica, por lo que este enfoque debe considerarse experimental. Las estadísticas de ajuste aplicadas a los historiales de encuentros agregados ofrecen un enfoque alternativo (MacKenzie y Bailey 2004), pero son difíciles de implementar cuando J * T es alto y hay valores faltantes o covariables continuas.

 chisq <- function(fm) { 
   umf <- getData(fm) 
   y <- getY(umf) 
   sr <- fm@sitesRemoved 
   if(length(sr)>0) 
     y <- y[-sr,,drop=FALSE] 
   fv <- fitted(fm, na.rm=TRUE) 
   y[is.na(fv)] <- NA 
   sum((y-fv)^2/(fv*(1-fv))) 
 }

set.seed(344) 
pb.gof <- parboot(m0, statistic=chisq, nsim=100)
plot(pb.gof)

La figura indica que, como se esperaba, el modelo de parámetro constante no se ajusta bien a los datos.

Modelo de Ocupación Dinámico con datos del mundo real

Este analisis es un tutorial del modelo de ocupación dinámico y usa parte de los datos y el analisis que forman parte del articulo:

Ahumada JA, Hurtado J, Lizcano D. 2013. Monitoring the Status and Trends of Tropical Forest Terrestrial Vertebrate Communities from Camera Trap Data: A Tool for Conservation. PLoS ONE. 8:e73707. doi:DOI: 10.1371/journal.pone.0073707.

Utilizamos datos de trampa de cámara que la Red de Tropical Ecology Assessment and Monitoring (TEAM) recolectados de forma regular a lo largo de un transecto del Volcán Barva en Costa Rica.

Mostramos cómo estos datos pueden usarse para calcular los indicadores temporales de las especies de mamíferos de interés en el área.

Se encontraron descensos en la ocupación de algunas especies.

Para este tutorial nos enfocaremos en Dasiprocta punctata

Información de la sesión en R.

print(sessionInfo(), locale = FALSE)
## R version 4.0.3 (2020-10-10)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 14393)
## 
## Matrix products: default
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] unmarked_1.1.1  lattice_0.20-44 ggplot2_3.3.6   markdown_1.1   
## 
## loaded via a namespace (and not attached):
##  [1] tidyselect_1.1.2 terra_1.5-21     xfun_0.31        bslib_0.4.0     
##  [5] purrr_0.3.4      splines_4.0.3    colorspace_2.0-3 vctrs_0.4.1     
##  [9] generics_0.1.3   htmltools_0.5.3  yaml_2.3.5       utf8_1.2.2      
## [13] rlang_1.0.4      nloptr_2.0.1     jquerylib_0.1.4  pillar_1.8.0    
## [17] glue_1.6.2       withr_2.5.0      DBI_1.1.2        sp_1.4-5        
## [21] lifecycle_1.0.1  plyr_1.8.7       stringr_1.4.0    munsell_0.5.0   
## [25] gtable_0.3.0     raster_3.5-15    codetools_0.2-18 evaluate_0.15   
## [29] labeling_0.4.2   knitr_1.39       fastmap_1.1.0    parallel_4.0.3  
## [33] fansi_1.0.3      highr_0.9        Rcpp_1.0.8.3     scales_1.2.0    
## [37] cachem_1.0.6     jsonlite_1.8.0   lme4_1.1-29      farver_2.1.0    
## [41] digest_0.6.29    stringi_1.7.8    dplyr_1.0.9      grid_4.0.3      
## [45] cli_3.3.0        tools_4.0.3      magrittr_2.0.3   sass_0.4.2      
## [49] tibble_3.1.8     pkgconfig_2.0.3  Matrix_1.3-4     MASS_7.3-54     
## [53] minqa_1.2.4      assertthat_0.2.1 rmarkdown_2.14   rstudioapi_0.13 
## [57] boot_1.3-28      R6_2.5.1         nlme_3.1-152     compiler_4.0.3