Taller SECEM

Presentación

En este espacio se encuentran los contenidos teóricos y prácticos para el taller Introducción a la modelización jerárquica: detectabilidad imperfecta en modelos de ocupación y N-mixture impartido en el XVI Congreso de la SECEM (6-9 de diciembre de 2023 en Granollers). En los siguientes enlaces se pueden descargar los materiales se utilizarán durante el taller:

• Presentaciones de los contenidos teóricos.
• Script en R y datos necesarios para las sesiones prácticas.

Se recomienda guardar todos estos archivos en una misma carpeta. A continuación se presenta de manera resumida la teoría del taller y posteriormente los casos prácticos.

Proceso ecológico y observacional

Cuando estudiamos la distribución o la abundancia (variables de estado) de las especies podemos distinguir entre

• Proceso ecológico: Mecanismo que da lugar a la variable de estado (distribución o abundancia en nuestro caso) mediante el cual un conjunto de covariables abióticas (temperatura, composición del suelo, etc.), bióticas (interacciones entre especies, etc.) o demográficas (tasas de reproducción, mortalidad, dispersión, etc.) originan que los organismos se distribuyan de una determinada forma en nuestro área de estudio.
• Proceso observacional: Mecanismo mediante el cual, a partir de una distribución o abundancia dada, emergen nuestros datos (observaciones) recogidos en el campo: fotos de una cámara de fototrampeo, conteos directos, registros de bioacústica, registro de indicios, muestreos de ADN ambiental, registros de ciencia ciudadana, etc.

Si no se tiene en cuenta ese proceso observacional corremos el riesgo de asumir que nuestros datos representan fielmente la realidad, lo cual podría hacernos obtener estimas de la distribución o la abundancia sesgadas y/o poco precisas.

Detectabilidad imperfecta

Uno de los sesgos más comunes en los procesos observacionales a la hora de estudiar la distribución y la abundancia de las especies es la detectabilidad imperfecta, esto es, la probabilidad de que una especie o un individuo no sea detectado estando presente en un lugar dado. Este fenómeno es muy habitual, ya que hay muchos factores que pueden afectar a la detectabilidad: el carácter elusivo de las especies, su fenología, los medios de los que disponemos para detectarlas, las condiciones metereológicas en el momento del muestreo, la experiencia del muestreador, el propio hábitat donde se encuentra la especie, etc. La detectabilidad imperfecta es especialmente importante cuando varía a través de nuestros sitios de muestreo, ya que podemos confundir el efecto de las covariables que afectan a la detectabilidad con las que afectan a la ecología de nuestra especie.

Modelos de ocupación

Los modelos de ocupación nos permiten estimar la probabilidad de presencia o distribución de una especie (esta será nuestra variable de estado) en nuestro área de estudio teniendo en cuenta la detectabilidad imperfecta a partir de muestreos repetidos registrando detecciones/no detecciones (observaciones). Para ello, asumiremos que el proceso de que una especie esté presente en el sitio \(i\) (\(z_{i} = 1\)) sigue una distribución de \(Bernoulli\) con probabilidad de éxito (probabilidad de ocupación) \(\psi_{i}\):

\[\begin{equation} \tag{Eq 1}\label{eq1} z_{i} \sim Bernoulli(\psi_{i}) \end{equation}\]

Esta probabilidad de ocupación a su vez puede ser modelizada por medio de covariables predictoras utilizando la maquinaria habitual de los modelos lineales generalizados, utilizando por ejemplo la función vínculo \(logit\) para asegurarnos de que el parámetro \(\psi_{i}\) toma valores entre 0 y 1, ya que se trata de una probabilidad (otras funciones vínculo como \(probit\) o \(cloglog\) son igualmente posibles):

\[\begin{equation} \tag{Eq 2}\label{eq2} logit(\psi_{i}) = \boldsymbol{\beta X_{i}} \end{equation}\]

donde \(\boldsymbol{\beta}\) representa el vector de que forman el intercepto \(\beta_{0}\) y los coeficientes o estimates de cada una de las covariables predictoras (\(\beta_{1}, \beta_{2}, ... \beta_{n}\)), mientras que \(\boldsymbol{X_{i}}\) representa el valor de cada una de las covariables predictoras (\(X1_{i}, X2_{i}, ... Xn_{i}\)) en el sitio \(i\).

Hasta aquí, el procedimiento es exáctamente idéntico a una regresión logística utilizada normalmente para modelizar la presencia/ausencia de las especies. La diferencia en los modelos de ocupación radica en que se añade un sub-modelo que describe el proceso observacional y se sirve del histórico de detecciones/no detecciones para estimar la probabilidad de detección de la especie en cada sitio \(i\):

\[\begin{equation} \tag{Eq 3}\label{eq3} h_{ij}|z_{i} \sim Bernoulli(z_{i}p_{ij}) \end{equation}\]

donde \(h_{ij}|z_{i}\) es el histórico de \(j\) detecciones/no detecciones dado que la especie esté realmente presente o ausente (1 / 0) en el sitio \(i\) (\(z_{i}\)), mientras que \(p_{ij}\) es la probabilidad de detección de la especie en el sitio \(i\) durante el muestreo \(j\). Una de las ventajas de los modelos de ocupación es que nos permiten a la vez modelizar esa probabilidad de detección \(p_{ij}\) por medio de covariables predictoras de la misma forma que hicimos con la ocupación:

\[\begin{equation} \tag{Eq 4}\label{eq4} logit(p_{ij}) = \boldsymbol{\beta X_{ij}} \end{equation}\]

Cabe destacar que en el proceso observacional pueden existir dos tipos de covariables predictoras: covariables que afecten al sitio \(i\), o covariables que afecten a la ocasión de muestreo \(j\). Puede ser que la detectabilidad de la especie se vea afectada por las características físicas de una determinada localidad \(i\) si por ejemplo nos encontramos en un terreno muy escarpado con mucho desnivel y esto nos impide localizar fácilmente a nuestra especie en comparación con localidades más planas. Pero también puede suceder que las condiciones meteorológicas en una determinada ocasión de muestreo impidan la detección de la especie (por ejemplo, momentos de mucha niebla) en comparación a otras ocasiones de muestreo en la misma localidad (días despejados). Nótese que esta es una característica propia del proceso observacional, ya que en los modelos de ocupación se asume que la variable de estado (presencia o ausencia real de la especie en una localidad) no varía durante la realización del estudio. Esta situación es igualmente aplicable a los modelos N-mixture que veremos a continuación.

De esta forma se puede apreciar la estructura jerárquica de los modelos de ocupación, ya que el proceso observacional (\ref{eq3} y \ref{eq4}) se construye sobre la base de que la especie esté realmente presente o ausente en cada sitio (proceso ecológico, \ref{eq1} y \ref{eq2}), pudiendo estar afectado independientemente por otras covariables.

Modelos N-mixture

Los modelos N-mixture son la versión para trabajar con conteos y abundancias de los modelos de ocupación. Nos permiten estimar la abundancia de una especie (esta será nuestra variable de estado) en nuestro área de estudio teniendo en cuenta la detectabilidad imperfecta a partir de muestreos repetidos registrando conteos de los individuos de nuestra especie (observaciones). Para ello, asumiremos que la abundancia de una especie en el sitio \(i\) (\(N_{i}\)) sigue una distribución de \(Poisson\) con intensidad o abundancia media \(\lambda_{i}\):

\[\begin{equation} \tag{Eq 5}\label{eq5} N_{i} \sim Poisson(\lambda_{i}) \end{equation}\]

De forma similar a los modelos de ocupación, esta intensidad o abundancia media \(\lambda_{i}\) puede modelizarse a través de covariables predictoras utilizando la maquinaria habitual en los modelos lineales generalizados utilizando esta vez la función vínculo \(log\) para que el valor de \(\lambda_{i}\) siempre sea positivo:

\[\begin{equation} \tag{Eq 6}\label{eq6} log(\lambda_{i}) = \boldsymbol{\beta X_{i}} \end{equation}\]

donde, análogamente al modelo de ocupación, \(\boldsymbol{\beta}\) representa el vector de que forman el intercepto \(\beta_{0}\) y los coeficientes o estimates de cada una de las covariables predictoras (\(\beta_{1}, \beta_{2}, ... \beta_{n}\)), mientras que \(\boldsymbol{X_{i}}\) representa el valor de cada una de las covariables predictoras (\(X1_{i}, X2_{i}, ... Xn_{i}\)) en el sitio \(i\).

Una vez más, hasta aquí el procedimiento sería idéntico a una regresión de Poisson, o un modelo lineal generalizado utilizando la familia de \(Poisson\). Sin embargo, ahora podemos añadir un nuevo sub-modelo que se encargue de modelizar el proceso observacional y tener en cuenta la detectabilidad imperfecta a la hora de contar los indivduos de nuestra especie. Para ello utilizaremos la distribución \(Binomial\) en la que cada uno de los individuos realmente presentes en cada sitio de muestreo \(i\) será un ensayo de “éxito o fracaso” que dependerá de una probabilidad de detección individual (\(p_{ij}\)):

\[\begin{equation} \tag{Eq 7}\label{eq7} C_{ij}|N_{i} \sim Binomial(N_{i},p_{ij}) \end{equation}\]

En este caso, \(C_{ij}\) será el resultado del conteo de individuos en el sitio \(i\) durante el muestreo \(j\), mientras que \(N_{i}\) será el número de individuos realmente presentes en el sitio \(i\). De igual forma que en los modelos de ocupación, la probabilidad de detección (esta vez a nivel de individuo), puede ser a su vez modelizada mediante covariables predictoras (nótese que volvemos a usar la función vínculo \(logit\) ya que la probabilidad debe encontrarse entre 0 y 1):

\[\begin{equation} \tag{Eq 8}\label{eq8} logit(p_{ij}) = \boldsymbol{\beta X_{ij}} \end{equation}\]

Entre las numerosas extensiones disponibles de este modelo general, una de las más básicas es la de asumir que la abundancia de nuestra especie sigue otra distribución que no sea la de \(Poisson\) determinada en la \ref{eq7}, como la \(\text{Binomial negativa}\) o la \(\text{Poisson inflada de ceros}\)

Práctica 1: Modelos de ocupación

Queremos modelizar la ocupación (probabilidad de presencia) de una especie a partir de datos de cámaras trampa. Nuestro área de estudio es una zona montañosa dividida en una malla de 50x50 celdas. Sabemos que nuestra especie utiliza mayormente el hábitat de matorral de los fondos de valle de nuestro área de estudio. Disponemos de un total de 100 cámaras de fototrampeo que hemos dispuesto aleatoriamente en nuestro área de estudio que han estado activas durante 7 noches utilizando una lata de sardinas como atrayente. De las 100 cámaras, 47 de ellas son de la marca A y 53 de la marca B.

# Instalamos y cargamos las librerías que vamos a utilizar.
#install.packages("terra")
#install.packages("unmarked")
#install.packages("AICcmodavg")
library(terra)
library(unmarked)
library(AICcmodavg)

# Leemos nuestros datos en formato CVS
datos <- read.csv("https://raw.githubusercontent.com/jabiologo/web/master/tutorials/tallerSECEM_files/occuDatos.csv")
# Si hemos descargado previamente los datos, podemos indicar la ruta (la carpeta
# en nuestro ordenador donde se localizan los archivos) y cargarlos desde ahí.
# datos <- read.csv("mi/ruta/occuDatos.csv")

# Le pedimos a R que nos muestre las primeras filas de nuestros datos.
head(datos)

##     id marca o1 o2 o3 o4 o5 o6 o7 dia1 dia2 dia3 dia4 dia5 dia6 dia7
## 1  593     A  0  0  0  0  0  0  0    1    2    3    4    5    6    7
## 2 1942     B  0  0  0  0  0  0  0    1    2    3    4    5    6    7
## 3 1945     A  0  0  0  1  0  0  0    1    2    3    4    5    6    7
## 4  253     B  0  0  0  0  0  0  0    1    2    3    4    5    6    7
## 5 2152     A  0  0  0  0  0  0  0    1    2    3    4    5    6    7
## 6 1528     A  0  0  0  0  0  0  0    1    2    3    4    5    6    7

# Cargamos las capas raster correspondientes a nuestras covariables predictoras,
# en este caso elevación y porcentaje de vegetación arbustiva.
elev <- rast("https://github.com/jabiologo/web/raw/master/tutorials/tallerSECEM_files/elev.tif")
arbu <- rast("https://github.com/jabiologo/web/raw/master/tutorials/tallerSECEM_files/arbu.tif")
# Si hemos descargado previamente las capas, podemos indicar la ruta (la carpeta
# en nuestro ordenador donde se localizan los archivos) y cargarlos desde ahí.
# elev <- rast("mi/ruta/elev.tif")
# arbu <- rast("mi/ruta/arbu.tif")

names(elev) <- "elev"
names(arbu) <- "arbu"

# Vamos a graficar estas capas para hacernos una idea de cómo son. Además,
# colocaremos encima las localizaciones de nuestros dos modelos de cámaras.
par(mfrow = c(1,2))
plot(elev, main = "Elevación")
points(xyFromCell(elev, datos$id[datos$marca == "A"]), col = "darkred", pch = 19)
points(xyFromCell(elev, datos$id[datos$marca == "B"]), col = "darkblue", pch = 19)
plot(arbu, main = "Porcentaje arbusto")

# En este paso estandarizaremos nuestras variables para que su media sea 0 y su
# desviación estándar 1. Esta es una práctica habitual en modelización que ayuda
# a ajustar los modelos de forma más efectiva. La fórmula es:
# valor escalado = (valor de la cov - valor medio de la cov)/ SD de la cov
datos$elev <- scale(elev[datos$id])[1:100]
datos$arbu <- scale(arbu[datos$id])[1:100]

# Volvemos a inspeccionar las primeras filas de nuestro juego de datos.
head(datos)

##     id marca o1 o2 o3 o4 o5 o6 o7 dia1 dia2 dia3 dia4 dia5 dia6 dia7
## 1  593     A  0  0  0  0  0  0  0    1    2    3    4    5    6    7
## 2 1942     B  0  0  0  0  0  0  0    1    2    3    4    5    6    7
## 3 1945     A  0  0  0  1  0  0  0    1    2    3    4    5    6    7
## 4  253     B  0  0  0  0  0  0  0    1    2    3    4    5    6    7
## 5 2152     A  0  0  0  0  0  0  0    1    2    3    4    5    6    7
## 6 1528     A  0  0  0  0  0  0  0    1    2    3    4    5    6    7
##          elev       arbu
## 1  0.55262104  1.4229936
## 2  0.37677981 -0.4942135
## 3  0.03575438  0.7839246
## 4 -0.80082361  0.9685444
## 5  0.29685197  0.3578785
## 6  0.35546572 -1.5735298

# A continuación prepararemos los datos para que sean entendidos por el paquete
# unmarked. Para ello debemos colocar de una forma determinada nuestras
# observaciones (detecciones/no detecciones en este caso), nuestras covariables
# predictoras. Recordad que las covariables predictoras que pueden estar
# relacionadas con los sitios de muestreos (siteCovs) o con las ocasiones de
# muestreo (obsCovs).
obserCov <- list(dia = as.matrix(datos[,10:16]))
datosUM <- unmarkedFrameOccu(as.matrix(datos[,3:9]), siteCovs=datos[,c(2,17,18)], obsCovs=obserCov)

## Warning: siteCovs contains characters. Converting them to factors.

# A continuación ajustaremos el modelo de ocupación utilizando la marca de la 
# cámara trampa y los días desde que se colocó la cámara y el cebo como 
# covariables que afectan a la detectabilidad (proceso observacional), y 
# elevación y porcentaje de cobertura arbustiva como covariables que afectan a 
# la ocupación o probabilidad de presencia (proceso ecológico).
m1 <- occu(~ marca  + dia ~ elev + arbu, data = datosUM)

#Inspeccionamos los resultados del modelo.
summary(m1)

## 
## Call:
## occu(formula = ~marca + dia ~ elev + arbu, data = datosUM)
## 
## Occupancy (logit-scale):
##             Estimate    SE      z  P(>|z|)
## (Intercept)  -0.6183 0.273 -2.266 0.023439
## elev         -1.0045 0.300 -3.354 0.000797
## arbu         -0.0257 0.254 -0.101 0.919423
## 
## Detection (logit-scale):
##             Estimate     SE      z  P(>|z|)
## (Intercept)    0.218 0.4549  0.478 6.32e-01
## marcaB         0.690 0.4026  1.713 8.67e-02
## dia           -0.443 0.0831 -5.331 9.79e-08
## 
## AIC: 368.7241 
## Number of sites: 100
## optim convergence code: 0
## optim iterations: 32 
## Bootstrap iterations: 0

# Inspeccionamos los efectos de la covariable de la cámara  y del procentaje 
# de arbusto.
plotEffects(m1, type = "det", covariate="marca")

plotEffects(m1, type = "state", covariate="arbu")

# También podemos predecir de forma espacialmente explícita nuestro modelo 
# de ocupación para inspeccionar la probabilidad de presencia de nuestra especie
# sobre un mapa. Para ello primero estandarizaremos los rásters de la misma
# forma que lo hicimos con nuestras variables.
elevScaled <- (elev - 958.29) / 187.6693
arbuScaled <- (arbu - 49.78) / 7.041493
dataPred <- c(elevScaled, arbuScaled)
names(dataPred) <- c("elev", "arbu")

# Después, podemos utilizar la función predict con nuestro modelo y las nuevas
# capas raster estandarizadas y graficar los resultados.
pred <- predict(m1, dataPred, type = "state")
plot(pred)

# Por último podemos explorar la bondad de ajuste de nuestro modelo, esto es,
# cuanto de bien predice nuestras observaciones el modelo que hemos construido.
mb.gof.test(m1, nsim = 100)

## 
## MacKenzie and Bailey goodness-of-fit for single-season occupancy model
## 
## Pearson chi-square table:
## 
##         Cohort Observed Expected Chi-square
## 0000000      0       67    66.67       0.00
## 0001000      0        5     1.05      14.94
## 0010000      0        1     1.63       0.24
## 0010100      0        1     0.31       1.54
## 0011010      0        1     0.07      13.26
## 0100000      0        1     2.54       0.93
## 0110000      0        1     1.17       0.02
## 0111110      0        1     0.01      67.15
## 1000000      0        6     3.95       1.06
## 1000001      0        2     0.31       9.24
## 1001000      0        1     1.17       0.02
## 1010000      0        1     1.82       0.37
## 1010100      0        1     0.39       0.97
## 1011010      0        1     0.09       9.15
## 1100000      0        2     2.83       0.24
## 1100010      0        1     0.39       0.97
## 1100100      0        1     0.60       0.26
## 1110000      0        2     1.46       0.20
## 1110010      0        1     0.22       2.78
## 1110100      0        1     0.34       1.27
## 1110110      0        1     0.05      16.36
## 1111000      0        1     0.53       0.41
## 
## Chi-square statistic = 153.8116 
## Number of bootstrap samples = 100
## P-value = 0.16
## 
## Quantiles of bootstrapped statistics:
##   0%  25%  50%  75% 100% 
##   51   79   95  131  273 
## 
## Estimate of c-hat = 1.4

¿Qué podemos decir de este modelo? En primer lugar, puede sorprendernos el coeficiente estimado para la covariable porcentaje de arbusto. Este coeficiente es ligeramente negativo (-0.02), pero es poco significativo, su p-value es mayor de 0.05 y el intervalo de confianza de este valor utilizando el error estándar incluye el valor 0, por lo que esta covariable predictora parece que tiene poco efecto en nuestra especie. Sin embargo, como investigadores estábamos seguros de que esta especie utilizaba el hábitat de matorral del fondo de valle de nuestra área de estudio…
Cuando graficamos los efectos de las covariables en nuestro modelo también podemos apreciar que la marca de cámara trampa B es bastante más efectiva que la A a la hora de detectar nuestra especie, y confirmamos que el porcentaje de arbusto apenas tiene efecto sobre la probabilidad de presencia de nuestra especie.
Por último podemos comprobar que la bondad de ajuste de nuestro modelo no es demasiado buena… el p-valor no es demasiado significativo pero atendiendo al valor de ĉ apreciamos algo de sobredispersión.
Quizás podríamos construir otro modelo modificando las covariables predictoras en los diferentes sub-modelos para los procesos ecológico y observacional y compararlo con el modelo 1 mediante su valor de AIC.

Práctica 2: Modelos N-mixture

En este caso queremos modelizar la abundancia (número de individuos) de una especie a partir de conteos realizados en muestreos repetidos. Nuestro área de estudio es una zona de transición agrícola-forestal en una malla de 50x50 celdas. Sabemos que nuestra especie habita principalmente en bosques y tiene querencia por temperaturas cálidas. Hemos muestreado un total de 100 sitios diferentes en los que hemos realizado conteos directos en tres ocasiones distintas. Como sabemos que nuestra especie es sensible a la temperatura, en cada muestreo tomamos nota de la temperatura que había en ese momento en ese lugar, y le dimos un valor categórico (temperatura baja, media o alta).

# Si hemos cargado los paquetes necesarios anteriormente, no es necesario volver
# a cargarlos. Los paquetes que utilizaremos terra, unmarked y AICcmodavg
# Leemos nuestros datos en formato CVS.
datos <- read.csv("https://github.com/jabiologo/web/raw/master/tutorials/tallerSECEM_files/nmixDatos.csv")
# Si hemos descargado previamente los datos, podemos indicar la ruta (la carpeta
# en nuestro ordenador donde se localizan los archivos) y cargarlos desde ahí
# datos <- read.csv("mi/ruta/nmixDatos.csv").

# Visualizamos las primeras filas de nuestros datos.
head(datos)

##     id o1 o2 o3    t1    t2    t3
## 1  553  0  0  0 media media media
## 2 1897  8  7  1  alta  alta  baja
## 3 1947  0  0  0 media  alta media
## 4 2283  3  1  3 media media  alta
## 5 1549 10 10 10 media media  alta
## 6  933  1  0  4  baja  baja media

# Cargamos las capas raster correspondientes a nuestras covariables predictoras,
# en este caso tipo de cobertura de vegetación y temperatura media de cada sitio.
cober <- rast("https://github.com/jabiologo/web/raw/master/tutorials/tallerSECEM_files/cober.tif")
tempe <- rast("https://github.com/jabiologo/web/raw/master/tutorials/tallerSECEM_files/tempe.tif")
# Si hemos descargado previamente las capas, podemos indicar la ruta (la carpeta
# en nuestro ordenador donde se localizan los archivos) y cargarlos desde ahí.
# cober <- rast("mi/ruta/cober.tif")
# tempe <- rast("mi/ruta/tempe.tif")
names(cober) <- "cober"
names(tempe) <- "tempe"

# Vamos a graficar estas capas para hacernos una idea de cómo son. Colocaremos 
# encima las localizaciones de nuestros conteos repetidos. Nótese que la 
# variable cobertura del suelo es categórica y cada númer corresponde a un tipo
# de cobertura:
# 1 = urbano
# 2 = agrícola
# 3 = transición
# 4 = bosque
par(mfrow = c(1,2))
plot(tempe, main = "Temperatura")
points(xyFromCell(tempe, datos$id), pch = 19, cex = 0.5)
plot(cober, main = "Cobertura del suelo")

# Estandarizamos la variable temperatura y convertimos a factor la variable de
# cobertura del suelo, ya que se trata de una variable categórica.
datos$tempe <- scale(tempe[datos$id])[1:100]
datos$cober <- as.factor(cober[datos$id][,1])

# Volvemos a visualizar las primeras filas de nuestro juego de datos.
head(datos)

##     id o1 o2 o3    t1    t2    t3      tempe cober
## 1  553  0  0  0 media media media -0.4727608     2
## 2 1897  8  7  1  alta  alta  baja -2.0753400     4
## 3 1947  0  0  0 media  alta media -1.4743728     3
## 4 2283  3  1  3 media media  alta -1.4743728     3
## 5 1549 10 10 10 media media  alta -1.4743728     4
## 6  933  1  0  4  baja  baja media  0.5288511     3

# Como hicimos anteriormente, vamos a organizar nuestros datos de una manera
# que entienda el paquete unmarked.
obserCov <- list(tOc = datos[,5:7])
datosUM <- unmarkedFramePCount(y = as.matrix(datos[,2:4]), siteCovs=datos[,8:9], obsCovs=obserCov)

## Warning: obsCovs contains characters. Converting them to factors.

# A continuación ajustaremos un modelo N-mixture utilizando la temperatura de la
# ocasión de muestreo y la temperatura media del sitio como covariables que
# afectan a la detectabilidad y la temperatura media del sitio y el tipo de
# cobertura del suelo como variables que afectan a la abundancia de la especie.
m1 <- pcount(~ tOc + tempe ~ tempe + cober, data = datosUM, K = 200)

# Inspeccionamos los resultados del modelo
summary(m1)

## 
## Call:
## pcount(formula = ~tOc + tempe ~ tempe + cober, data = datosUM, 
##     K = 200)
## 
## Abundance (log-scale):
##             Estimate     SE     z  P(>|z|)
## (Intercept)    0.484 0.5761  0.84 4.01e-01
## tempe          1.079 0.0718 15.02 5.33e-51
## cober2        -2.115 0.6444 -3.28 1.03e-03
## cober3         0.887 0.5810  1.53 1.27e-01
## cober4         3.816 0.5836  6.54 6.17e-11
## 
## Detection (logit-scale):
##             Estimate    SE       z  P(>|z|)
## (Intercept)    1.313 0.281   4.666 3.07e-06
## tOcbaja       -3.428 0.262 -13.102 3.19e-39
## tOcmedia      -0.550 0.196  -2.815 4.88e-03
## tempe         -0.059 0.139  -0.424 6.72e-01
## 
## AIC: 590.2252 
## Number of sites: 100
## optim convergence code: 0
## optim iterations: 86 
## Bootstrap iterations: 0

# Inspeccionamos algunos efectos de algunas covariables
plotEffects(m1, type = "det", covariate="tempe")
plotEffects(m1, type = "state", covariate="tempe")

# Preparamos las variables para poder realizar las predicciones
tempeScaled <- (tempe[] - 26.22) / 5.02
dataPred <- data.frame(cbind(tempeScaled, cober[]))
names(dataPred) <- c("tempe", "cober")
dataPred$cober <- as.factor(dataPred$cober)

# Realizamos las predicciones del modelo y las graficamos
pred <- predict(m1, dataPred, type = "state")
predRas <- rast(nrows=50, ncols=50, xmin=0, xmax=50, ymin = 0, ymax = 50)
predRas[] <- pred$Predicted
plot(predRas)

# Por último, podemos analizar la bondad de ajuste de nuestro modelo.
Nmix.gof.test(m1, nsim = 100)

## 
## Chi-square goodness-of-fit for N-mixture model of 'unmarkedFitPCount' class
## 
## Observed chi-square statistic = 207.3841 
## Number of bootstrap samples = 100
## P-value = 0.93
## 
## Quantiles of bootstrapped statistics:
##   0%  25%  50%  75% 100% 
##  165  233  261  305  949 
## 
## Estimate of c-hat = 0.72

Al igual que en el caso anterior, nos podemos preguntar muchas cosas sobre este modelo. ¿Qué efecto parece que tiene la temperatura media del sitio de muestreo sobre la detectabilidad de los individuos de nuestra especie? ¿Cómo podríamos mejorar este modelo?
Hay otras preguntas sobre las que podemos reflexionar comparando los modelos de ocupación y los modelos N-mixture:

• ¿Las detectabilidades \(p_{ij}\) que estamos calculando en cada modelo son la misma métrica?
• ¿Afecta de igual forma la asunción de población cerrada en los modelos de ocupación y en los modelos N-mixture?
  

Extensiones

A continuación se mencionan algunas extensiones a estos modelos y otras lecturas recomendadas:

• Modelos de ocupación con más de un estado, interesantes para modelizar diferentes categorías de ocupación, como ausencia/presencia/presencia con crías o ausencia/presencia/presencia con infección. Función occuMS del paquete unmarked (J. Royle et al. 2007).
• Modelos de ocupación teniendo en cuenta falsos positivos, interesantes cuando consideramos que la probabilidad de cometer falsos positivos es alta, como a la hora de identificar una especie a partir de contactos indirectos o genéticamente cuando la incertidumbre es alta (J. A. Royle and Link 2006). Función occuFPdel paquete unmarked.
• Modelos de Royle-Nichols, para estimas de abundancias a partir de detecciones/no detecciones de una especie. Función occuRN del paquete unmarked (J. A. Royle and Nichols 2003). Un ejemplo en Fernández-López et al. (2022).
• Modelos dinámicos, interesantes para dinámicas poblacionales temporales o cuando tenemos certeza de que se viola la asunción de población cerrada. Funciones colext() (MacKenzie et al. 2003) y pcountOpen() (Dail and Madsen 2011) del paquete unmarked.
• Modelos de ocupación multi específicos. Función occuMulti del paquete unmarked.
• Modelos integrados verosimilitud conjunta, cuando queremos combinar más de una fuente de datos en un único modelo. Tutorial utilizando NIMBLE en Fernández-López, Acevedo, and Gimenez (2023).

Go Bayesian!

Entre otras ventajas, la Inferencia Bayesiana proporciona una mayor flexibilidad y permite ajustar modelos mucho más complejos, que en ocasiones son necesarios. A continuación se mencionan algunas herramientas que nos permiten ajustar modelos jerárqucos mediante Inferencia Bayesiana con algunas extensiones, como efectos aleatorios, modelos espaciales, etc.

• paquete de R hSDM
• paquete de R ubms
• paquetes de R spOccupancy y spAbundance
• INLA
• NIMBLE, JAGS, BUGS, Stan, etc.

Referencias

Dail, D, and L Madsen. 2011. “Models for Estimating Abundance from Repeated Counts of an Open Metapopulation.” Biometrics 67: 577–87. https://doi.org/https://doi.org/10.1111/j.1541-0420.2010.01465.x.

Fernández-López, J, P Acevedo, and O Gimenez. 2023. “La Unión Hace La Fuerza: Modelos de Distribución de Especies Integrando Diferentes Fuentes de Datos.” Ecosistemas 32 (1): 2527–27. https://doi.org/https://doi.org/10.7818/ECOS.2527.

Fernández-López, J, JA Blanco-Aguiar, J Vicente, and P Acevedo. 2022. “Can We Model Distribution of Population Abundance from Wildlife–Vehicles Collision Data?” Ecography e06113. https://doi.org/10.1111/ecog.06113.

MacKenzie, D I, J A 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. https://doi.org/https://doi.org/10.1890/02-3090.

Royle, J A, and W A Link. 2006. “N-Mixture Models for Estimating Population Size from Spatially Replicated Counts.” Ecology 87: 35–841. https://doi.org/https://doi.org/10.1111/j.0006-341X.2004.00142.x.

Royle, J A, and J D Nichols. 2003. “N-Mixture Models for Estimating Population Size from Spatially Replicated Counts.” Ecology 84(3): 777–90. https://doi.org/ https://doi.org/10.1890/0012-9658(2003)084[0777:EAFRPA]2.0.CO;2.

Royle, JD, J D Hines, D I Mackenzie, M E Seamans, and R J Guitierrez. 2007. “Occupancy Estimation and Modeling with Multiple States and State Uncertainty.” Ecology 88: 1395–1400. https://doi.org/https://doi.org/10.1890/06-1474.