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)
## Warning: package 'unmarked' was built under R version 3.6.3
## Loading required package: lattice
## Warning: package 'lattice' was built under R version 3.6.3
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.

Información de la sesión en R.

print(sessionInfo(), locale = FALSE)
## R version 3.6.1 (2019-07-05)
## 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.0.1  lattice_0.20-41
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_1.0.5       codetools_0.2-16 digest_0.6.25    MASS_7.3-51.6   
##  [5] grid_3.6.1       plyr_1.8.6       magrittr_1.5     evaluate_0.14   
##  [9] rlang_0.4.6      stringi_1.4.6    sp_1.4-2         raster_3.3-7    
## [13] rmarkdown_2.3    tools_3.6.1      stringr_1.4.0    parallel_3.6.1  
## [17] xfun_0.15        yaml_2.2.1       compiler_3.6.1   htmltools_0.5.0 
## [21] knitr_1.29