Modelo de ocupaciĆ³n: Willow Tit

Objetivos:

  1. Hacer una regresiĆ³n logĆ­stica Bayesiana para modelar la probabilidad de ocupaciĆ³n de una especie en funciĆ³n de covariables, utilizando JAGS, bajo dos escenarios:
  • Asumiendo que no hay errores de detecciĆ³n.
  • Asumiendo detecciĆ³n imperfecta, incorporando la probabilidad de detecciĆ³n en modelos de ocupaciĆ³n.

Vamos a trabajar con datos de presencia/ausencia de ā€œWillow Titā€ (o Carbonero Montano, Poecile montanus) en puntos de muestreo de Suiza (ver Royle y Dorazio 2008, capĆ­tulo 3). Los datos provienen de un muestreo anual que se realiza durante la temporada reproductiva de la especie de interĆ©s. Las unidades de muestreo son cuadrantes de \(1\) km2, en total se relevaron \(237\) cuadrantes, que fueron visitados entre \(2\) y \(3\) veces dentro de la temporada. En cada uno de estos cuadrantes se registrĆ³ si la especie era detectada por el observador, ya sea de manera visual o acĆŗstica (valor = \(1\)), o no era detectada (valor = \(0\)).

A continuaciĆ³n vamos a cargamos los datos y usar el comando str para ver de quĆ© se trata

datos <- read.table("https://github.com/jmmorales/cursoMyD/raw/master/Data/willow.txt", header = TRUE)

str(datos)
## 'data.frame':    237 obs. of  5 variables:
##  $ y.1   : int  0 0 0 0 0 0 0 0 1 0 ...
##  $ y.2   : int  0 0 0 0 0 0 0 0 1 0 ...
##  $ y.3   : int  0 0 0 0 0 0 0 0 0 1 ...
##  $ elev  : int  420 450 1050 1110 510 630 590 530 1140 770 ...
##  $ forest: int  3 21 32 35 2 60 5 13 50 57 ...

En la variable y.1 estĆ”n los datos de la primera visita, en y.2 la segunda visita y en y.3 la tercer visita. En cada uno de estos casos, si un sitio en particular no pudo ser visitado, se registra como NA, si se detectĆ³ a la especie 1, y si no se observĆ³ la especie en esa visita un 0. Por otro lado, elev es la altitud (en metros sobre el nivel del mar) de cada sitio y forest es el porcentaje de cobertura de bosque.

El objetivo es relacionar la presencia de esta especie de pĆ”jaro con la cobertura de bosque y la altitud. Centramos y estandarizamos las dos covariables continuas para que sean de magnitudes comparables y creamos una nueva covariable, (el cuadrado de la altitud) porque sospechamos que hay una altitud Ć³ptima para esta especie.

alti  <- as.vector(scale(datos$elev, center = TRUE))
bosq  <- as.vector(scale(datos$forest, center = TRUE))
alti2 <- alti * alti

# Juntamos los datos de las visitas en un objeto que vamos a llamar "ym"
ym <- as.matrix(cbind(datos$y.1, datos$y.2, datos$y.3))  

#Contamos cuantas veces visitamos cada uno de los sitios
J <- rowSums(!is.na(ym)) 

# Creamos un vector M con el nĆŗmero total de sitios visitados    
M <- nrow(ym)  

# Contamos el nĆŗmero de veces que encontramos a la especie en cada sitio, 
# para generar el vector de observaciones.
y <- rowSums(ym, na.rm = TRUE)

Modelo de DetecciĆ³n Perfecta:

Ahora vamos a hacer una regresiĆ³n logĆ­stica Bayesiana, asumiendo que no hay error de observaciĆ³n. Esto quiere decir que, cuando no encontramos a la especie en un sitio es porque no estaba presente (Esto se conoce como ā€œDetecciĆ³n perfectaā€). Cuando esto ocurre, nuestra observaciĆ³n (que llamamos y) es igual al verdadero proceso de ocupaciĆ³n de la especie (que vamos a llamar z). La probabilidad de ocurrencia, que vamos a llamar \(\psi\), es el valor esperado de la variable de ocupaciĆ³n.

cat(file = "WillowDP.bug",
"
model { 
  for(i in 1: M){      
    y[i] ~ dbin(psi[i], J[i]) 
    logit(psi[i]) <- b0 + b1 * alti[i] + b2 * alti2[i] + b3 * bosq[i] 
  }

  b0 ~ dnorm(0, 0.01) 
  b1 ~ dnorm(0, 0.01) 
  b2 ~ dnorm(0, 0.01) 
  b3 ~ dnorm(0, 0.01) 

  logit(psi0) <- b0
}
")

Como de costumbre, tenemos que armar una lista de los datos que le vamos a pasar al JAGS, escribir una funciĆ³n para los valores iniciales de las cadenas Markovianas, definir la lista de parĆ”metros que queremos guardar y los nĆŗmeros de iteraciones, ā€œthinningā€, burnin y cadenas.

data <- list ("y", "J", "M", "bosq", "alti", "alti2")

inits <- function(){
  list ( b0 = rnorm(1),
         b1 = rnorm(1),
         b2 = rnorm(1),
         b3 = rnorm(1))
}

parameters <- c("b0", "b1", "b2", "b3", "psi0") 

ni <- 5000  
nt <- 4     
nb <- 1000  
nc <- 3

Finalmente podemos llamar a JAGS

library(jagsUI)

wb.sim <- jags(data, inits, parameters, "willowDP.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: 237
##    Unobserved stochastic nodes: 4
##    Total graph size: 1999
## 
## 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.

Siempre que ajustamos un modelo usando MCMC tenemos que asegurarnos que las cadenas hayan convergido y que el tamaƱo efectivo de las muestras de las posteriores sean razonables:

print(wb.sim)
## JAGS output for model 'willowDP.bug', generated by jagsUI.
## Estimates based on 3 chains of 5000 iterations,
## adaptation = 100 iterations (sufficient),
## burn-in = 1000 iterations and thin rate = 4,
## yielding 3000 total samples from the joint posterior. 
## MCMC ran for 0.096 minutes at time 2021-12-08 14:44:03.
## 
##             mean    sd    2.5%     50%   97.5% overlap0 f  Rhat n.eff
## b0        -0.874 0.174  -1.216  -0.874  -0.527    FALSE 1 1.019   129
## b1         2.106 0.196   1.767   2.089   2.534    FALSE 1 1.107    30
## b2        -0.959 0.173  -1.331  -0.956  -0.630    FALSE 1 1.048    49
## b3         0.842 0.135   0.566   0.844   1.096    FALSE 1 1.013   261
## psi0       0.296 0.036   0.229   0.295   0.371    FALSE 1 1.021   122
## deviance 456.316 2.639 453.036 455.737 462.954    FALSE 1 1.012   193
## 
## **WARNING** Rhat values indicate convergence failure. 
## 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 = 3.4 and DIC = 459.764 
## DIC is an estimate of expected predictive error (lower is better).

Vemos que tanto la elevaciĆ³n como la cobertura de bosques tienen coeficientes positivos y distintos de cero mientras que la elevaciĆ³n al cuadrado tiene un coeficiente negativo (Implica que hay un mĆ”ximo de probabilidad de ocurrencia a elevaciones intermedias). Vemos tambiĆ©n que la probabilidad de ocurrencia para condiciones ambientales promedio es de 0.296.

Podemos graficar cĆ³mo cambian las probabilidades de ocurrencia predichas a partir de las posteriores de los coeficientes

op <- par(mfrow = c(1, 2), mar = c(5, 5, 3, 2) + 0.1, mgp = c(3.5, 1, 0), cex.lab = 1.3 , font.lab = 2, cex.axis = 1.3, bty = "n", las = 1)

plot(datos$elev, plogis(wb.sim$mean$b0 + wb.sim$mean$b1 * alti + wb.sim$mean$b2 * alti2), pch = 21, bg = "gray", xlab = "ElevaciĆ³n", ylab = "Probabilidad de Ocurrencia")

plot(datos$forest, plogis(wb.sim$mean$b0 + wb.sim$mean$b3 * bosq), pch = 21, bg = "gray", xlab = "Cobertura de bosque", ylab = "Probabilidad de Ocurrencia")

par(op)

Modelo de DetecciĆ³n Imperfecta:

Si consideramos que es posible que no se detecte a la especie en un sitio aunque estĆ© presente, podemos separar la presencia/ausencia ā€œobservadaā€ de la ā€œverdaderaā€, y distinguir entre la ā€œdetectabilidadā€ y la ā€œocurrenciaā€.

Cuando la detecciĆ³n es imperfecta, los valores \(0\) en nuestros datos pueden provenir de dos fuentes: (1) error de muestreo, o (2) ausencias reales. En la prĆ”ctica no podemos saber cuĆ”l es la fuente que originĆ³ un valor \(0\) en los datos, por lo que tenemos que diferenciar entre la observaciĆ³n y el verdadero proceso de ocupaciĆ³n de la especie, y explicitar la forma en la que se conectan ambos procesos.

Nuestra observaciĆ³n (y), depende de la probabilidad de detecciĆ³n (que llamaremos p). El verdadero proceso de ocupaciĆ³n de la especie (z), pasa a ser una variable oculta, que se expresa como una variable de distribuciĆ³n Bernoulli (una Binomial con \(N = 1\)) cuya probabilidad de Ć©xito (probabilidad de ocupaciĆ³n, \(\psi\)) depende de covariables como en el caso anterior.

La conexiĆ³n entre esa variable oculta z con las observaciones ocurre a travĆ©s de la probabilidad de detecciĆ³n p. Esto se hace definiendo a la observaciĆ³n (y) como un proceso Binomial, cuya probabilidad de Ć©xito es igual al producto entre z y p (ver como se expresa en el modelo BUGS de abajo).

Esta multiplicaciĆ³n permite que cuando z vale \(0\) (es decir, la especie estĆ” ausente), la probabilidad de Ć©xito del proceso de observaciĆ³n se hace \(0\) y tenemos un ā€œcero verdaderoā€ en los datos (cero generado por la ausencia real). Cuando z vale \(1\), entonces la probabilidad de Ć©xito de la binomial de las observaciones estĆ” dada por la probabilidad de detecciĆ³n p y podemos tener falsos negativos (ceros generados por la ā€œno detecciĆ³nā€ en el muestreo).

Esta posibilidad de separar el proceso ecolĆ³gico (presencia de una especie en funciĆ³n de variables ambientales en este caso) del proceso de observaciĆ³n, es una forma de modelado jerĆ”rquico o multinivel que resulta relativamente fĆ”cil de realizar con herramientas de anĆ”lisis Bayesiano. No es poca cosaā€¦

cat(file = "willow_DI.bug",
    "
 #Likelihood:
    model { 
     for(i in 1: M){ 
       y[i] ~ dbin(pex[i], J[i]) 
       pex[i] <- z[i] * p   
       z[i] ~ dbin(psi[i], 1)      
       logit(psi[i]) <- b0 + b1 * alti[i] + b2 * alti2[i] + b3 * bosq[i] 
    } 
  #Previas: 
    p  ~ dunif(0,1)                
    b0 ~ dnorm(0, 0.01) 
    b1 ~ dnorm(0, 0.01) 
    b2 ~ dnorm(0, 0.01) 
    b3 ~ dnorm(0, 0.01) 
    logit(psi0) <- b0
    }
    ")

Los datos que usamos para este modelo son los mismos que en el caso anterior. En la funciĆ³n inits tenemos que incluir valores iniciales para z y para p.

data <- list ( "y", "J", "M", "bosq", "alti", "alti2")
inits <- function(){
  list (z = rbinom(M, 1, 1),
        b0 = 4,
        p = runif(1),
        b1 = rnorm(1),
        b2 = rnorm(1),
        b3 = rnorm(1))
}

parameters <- c("p", "b0", "b1", "b2", "b3", "psi0")  
ni <- 5000  
nt <- 4     
nb <- 1000  
nc <- 3     

w.sim <- jags(data, inits, parameters, "willow_DI.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: 237
##    Unobserved stochastic nodes: 242
##    Total graph size: 2713
## 
## 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(w.sim)
## JAGS output for model 'willow_DI.bug', generated by jagsUI.
## Estimates based on 3 chains of 5000 iterations,
## adaptation = 100 iterations (sufficient),
## burn-in = 1000 iterations and thin rate = 4,
## yielding 3000 total samples from the joint posterior. 
## MCMC ran for 0.164 minutes at time 2021-12-08 14:44:10.
## 
##             mean     sd    2.5%     50%   97.5% overlap0     f  Rhat n.eff
## p          0.787  0.029   0.723   0.788   0.841    FALSE 1.000 1.001  1981
## b0        -0.207  0.274  -0.740  -0.208   0.345     TRUE 0.774 1.000  3000
## b1         2.065  0.303   1.495   2.054   2.707    FALSE 1.000 1.002   752
## b2        -1.147  0.277  -1.740  -1.141  -0.635    FALSE 1.000 1.000  3000
## b3         0.874  0.237   0.414   0.867   1.351    FALSE 1.000 1.000  2964
## psi0       0.449  0.067   0.323   0.448   0.585    FALSE 1.000 1.000  3000
## deviance 172.761 10.354 160.638 170.612 197.326    FALSE 1.000 1.002  3000
## 
## 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 = 53.6 and DIC = 226.373 
## DIC is an estimate of expected predictive error (lower is better).

Vemos que la probabilidad de detecciĆ³n es relativamente alta y que las variables ambientales tienen coeficientes similares a los de modelo anterior. La probabilidad de ocurrencia en condiciones promedio aumentĆ³.

Para graficar la posterior de la probabilidad de ocurrencia en condiciones promedio hacemos:

op <- par(cex.lab = 1.3 , font.lab = 2, cex.axis = 1.3, bty = "n", las = 1)
plot(density(w.sim$sims.list$psi0), main = "", xlab = "Pr. ocupaciĆ³n", lwd = 2)

par(op)

Podemos comparar grĆ”ficamente cĆ³mo los dos modelos producen distintos valores esperados de probabilidad de ocupaciĆ³n (el modelo de detecciĆ³n imperfecta estĆ” en rojo).

op <- par(mfrow = c(1, 2), mar = c(5, 5, 3, 2) + 0.1, mgp = c(3.5, 1, 0), cex.lab = 1.3 , font.lab = 2, cex.axis = 1.3, bty = "n", las = 1)

plot(datos$elev, plogis(wb.sim$mean$b0 + wb.sim$mean$b1 * alti + wb.sim$mean$b2 * alti2), pch = 21, bg = "gray", xlab = "ElevaciĆ³n", ylab = "Probabilidad de Ocurrencia", ylim = c(0, 1)) 

points(datos$elev, plogis(w.sim$mean$b0+w.sim$mean$b1*alti + w.sim$mean$b2*alti2), col = 2) 

plot(datos$forest, plogis(wb.sim$mean$b0+wb.sim$mean$b3*bosq), pch = 21, bg = "gray", xlab="Cobertura de bosque", ylab="Probabilidad de Ocurrencia", ylim = c(0, 1))
points(datos$forest, plogis(w.sim$mean$b0+w.sim$mean$b3*bosq), col = 2)

par(op)

Juan Manuel Morales. 6 de Septiembre del 2015. ƚltima actualizaciĆ³n: 2021-12-08