Movimiento de Elk (Cervus canadensis) Reintroducido en Bancroft, Ontario

Objetivos:

  1. Ajustar diferentes modelos de movimiento del tipo Random Walk (RW, Caminata aleatoria).
  2. Usar Posterior Predictive Checks para evaluar el ajuste de los modelos a los datos.

Vamos a ver ejemplos de un par de modelos de random walk (RW) ajustados a datos de trayectorias registradas con GPS. El trabajo completo estĆ” aquĆ­.

Vamos a cargar los datos y reordenarlos un poco.

elk <- read.table("https://sites.google.com/site/modelosydatos/elkGPS1.txt", 
    header = TRUE)

Esta base de datos contiene las coordenadas proyectadas (Easting: X, Northing: Y) para un Ćŗnico individuo (Id: elk-287). Cada una de estas coordenadas tiene asociados diferentes datos: (1) el tipo de hĆ”bitat al que corresponden las coordenadas, (2) la distancia entre dos coordenadas sucesivas, (3) los Ć”ngulos de giro, (4:13) son distancias desde el punto a diferentes componentes del paisaje, como la distancia al agua, por ejemplo.

Para ajustar los modelos de movimiento, en este caso las variables que nos interesan de la base de datos son los pasos (que vamos a llamar steps) y los Ć”ngulos de giro (que vamos a llamar turns). Los steps son las distancias (en km) entre posiciones consecutivas de los registros de GPS (uno para cada medianoche en este caso), y los turns son los Ć”ngulos de giro entre dos steps consecutivos. Tanto los steps como los turns se estiman a partir de dos datos de GPS consecutivos, el primer valor y el Ćŗltimo son NAs. Filtramos esos datos que son NA y definimos las dos variables:

step <- elk$km_day[2:163]  # Pasos
turn <- elk$turns[2:163]  # Angulos
n <- length(step)  # NĆŗmero de observaciones

Ahora vamos a ajustar una distribuciĆ³n Weibull a los steps y una Wrapped-Cauchy a los turns. Pero tenemos que considerar dos cosas: (1) Los ceros no estĆ”n definidos en la Weibull, asĆ­ que tenemos que reemplazarlos por valores pequeƱos. (2) La distribuciĆ³n Wrapped-Cauchy no estĆ” definida en JAGS, asĆ­ que tenemos que usar un truco para incluir esta distribuciĆ³n en el anĆ”lisis.

any(step == 0, na.rm = TRUE)
## [1] TRUE
step <- ifelse(step == 0, 0.01, step)

Ahora escribimos un modelo en lenguaje BUGS para un RW simple donde para cada movimiento obtenemos una distancia de desplazamiento de una Weibull y un valor de giro (respecto a la direcciĆ³n de movimiento del paso anterior) de una distribuciĆ³n Wrapped Cauchy. El objetivo es estimar los parĆ”metros de la Weibull y de la Wrapped Cauchy.

cat(file = "single.bug", "  model{
    for (t in 1:n) {
    # Likelihood para los steps
    step[t] ~ dweib(b, a)\t#b es el parƔmetro de forma
    
    # Likelihood para los turns
    ones[t] ~ dbern(wc[t])
    wc[t] <- (1/(2*Pi)*(1-pow(rho,2))/(1+pow(rho,2)-2*rho*cos(turn[t]-mu)))/100
    turn[t] ~ dunif(-3.14159265359, 3.14159265359)
  }
  
  # Previas
  a ~ dnorm(0, 0.01)T(0, ) 
  b ~ dnorm(0, 0.01)T(0, ) 
  
  rho ~ dunif(0, 1) 
  mu ~ dunif(-3.1416, 3.1416)
  
  Pi <- 3.14159265359
}")

Como siempre, antes de ajustar el modelo tenemos que definir quĆ© datos le vamos a pasar a JAGS, una funciĆ³n para generar valores iniciales de las cadenas Markovianas, y cuĆ”ntas iteraciones queremos correr

ones <- numeric(n) + 1
# Esta variable es para el 'truco' para definir la Wrapped-Cauchy.

data.single <- list("step", "turn", "n", "ones")

inits.single <- function() list(a = runif(1, 0.1, 2), b = runif(1, 0.4, 3), 
    mu = runif(1, 0, pi), rho = runif(1))

params.single <- c("a", "b", "mu", "rho")

ni <- 2000
nt <- 1
nb <- 1000
nc <- 3

Ahora usamos la funciĆ³n jags

library(jagsUI)

single.sim <- jags(data.single, inits.single, params.single, model.file = "single.bug", 
    n.chains = nc, n.burnin = nb, n.thin = nt, n.iter = ni)
## 
## Processing function input....... 
## 
## Done. 
##  
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 486
##    Unobserved stochastic nodes: 4
##    Total graph size: 2982
## 
## 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(single.sim)
## JAGS output for model 'single.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 0.068 minutes at time 2019-05-30 08:59:07.
## 
##              mean    sd     2.5%      50%    97.5% overlap0     f  Rhat
## a           1.166 0.093    0.995    1.162    1.354    FALSE 1.000 1.000
## b           0.577 0.032    0.512    0.577    0.642    FALSE 1.000 1.001
## mu          1.544 1.206   -2.631    1.828    2.934     TRUE 0.914 1.048
## rho         0.088 0.052    0.004    0.085    0.191    FALSE 1.000 1.010
## deviance 2986.988 2.717 2983.452 2986.488 2993.889    FALSE 1.000 1.002
##          n.eff
## a         3000
## b         1498
## mu          93
## rho        224
## deviance   974
## 
## 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 = 3.7 and DIC = 2990.673 
## DIC is an estimate of expected predictive error (lower is better).

Posterior Predictive Checking

Queremos ver si datos simulados con este modelo tienen caracterĆ­sticas similares a los datos observados.

Lo primero que vamos a hacer es simular pasos y Ɣngulos de giro, con las funciones rweibul y rwrpcauchy respectivamente, usando las distribuciones posteriores de los parƔmetros estimados con el modelo.

Algo que siempre hay que tener en cuanta cuando hacemos este tipo de simulaciones, es asegurarnos que R y BUGS usen la misma formulaciĆ³n para las distribuciones. En este caso, BUGS y R usan distintas parameterizaciones de la Weibull y entonces tenemos que hacer una transformaciĆ³n antes de simular valores de la Weibull.

R usa:

\[ f(x) = \frac{b}{a} \left( \frac{x}{a} \right)^{b - 1} \exp \left( - \left( \frac{x}{a} \right)^b \right) \]

Mientras que BUGS usa:

\[ f(x) = b a x^{ b - 1} \exp \left(- a x^{b} \right) \]

De manera que para ā€œtraducirā€ los valores de escala que estima BUGS a los que vamos a usar en R, tenemos que hacer \(\text{escala} = \left( \frac{1}{a} \right)^b\)

library(circular)  #Para usar la funciĆ³n rwrpcauchy()

n.sims <- single.sim$mcmc.info$n.samples  #Definimos el nĆŗmero de simulaciones
pps.single <- matrix(NA, n.sims, n)  # AcĆ” vamos a guardar los steps simulados
ppt.single <- matrix(NA, n.sims, n)  # AcĆ” vamos a guardar los turns simulados

a <- single.sim$sims.list$a
b <- single.sim$sims.list$b
mu <- single.sim$sims.list$mu
rho <- single.sim$sims.list$rho

# Simulamos los steps y los turns:
for (i in 1:n.sims) {
    pps.single[i, ] <- rweibull(n, shape = b[i], scale = (1/a[i])^b[i])
    ppt.single[i, ] <- rwrappedcauchy(n, mu = circular(mu[i]), rho = rho[i])
}

El segundo paso consiste en definir quĆ© mĆ©trica queremos usar para comparar los datos reales con los datos simulados a partir del modelo y las posteriores. En este ejemplo, nos interesa comparar visualmente la autocorrelaciĆ³n entre los steps observados y los simulados. Vamos a aplicar una funciĆ³n de autocorrelaciĆ³n (ACF por Auto-Correlation Function) a los datos simulados (lo vamos a llamar ppac.s), y a los datos observados (lo vamos a llamar oac).

nlags <- 61  #Definimos el nĆŗmero de lags que queremos usar en la ACF.
ppac.s <- matrix(NA, n.sims, nlags)  #AcĆ” vamos a guardar los valores de
# autocorrelaciĆ³n de los steps simulados.

# Estimamos la autocorrelaciĆ³n de los steps simulados:
for (i in 1:n.sims) {
    ppac.s[i, ] <- acf(log(pps.single[i, ]), lag.max = (nlags - 1), plot = FALSE)$acf
}

# Estimamos la autocorrelaciĆ³n de los steps observados
oac <- acf(log(step), lag.max = (nlags - 1), plot = FALSE)

library(coda)
hpd <- HPDinterval(as.mcmc(ppac.s))

# Graficamos para evaluar visualmente el ajuste del modelo a los datos:
op <- par(las = 1, cex.lab = 1.5, cex.axis = 1.3)
plot(oac$lag, oac$acf, type = "b", lwd = 2, col = 2, xlab = "lag", ylab = "autocorrelation")
lines(oac$lag, hpd[, 1], lwd = 2)
lines(oac$lag, hpd[, 2], lwd = 2)

par(op)

En color rojo vemos los valores de autocorrelaciĆ³n observados, y entre las lĆ­neas negras el intervalo de credibilidad de la autocorrelaciĆ³n en los datos simulados. Si el modelo ajustase bien a los datos, el intervalo de credibilidad deberĆ­a incluir los valores observados aproximadamente el \(95\)% de las veces. Vemos que el modelo de random walk simple no es capaz de producir un patrĆ³n de autocorrelaciĆ³n en distancia de desplazamiento como el observado.

Si graficamos la trayectoria del individuo, podemos ver que en algunas secciones de la trayectoria se compone de puntos que estƔn mƔs cercanos entre sƭ (steps mƔs cortos), mientras que en otras secciones los puntos estƔn mƔs espaciados (steps mƔs largos). En general, los steps mƔs cortos tienen giros mƔs amplios, y los steps mƔs largos estƔn asociados a giros mƔs pequeƱos.

plot(elk$Easting, elk$Northing, type = "o", asp = 1)

El prĆ³ximo paso es ver si mejora el ajuste a los datos con un modelo donde el animal alterna entre dos tipos de random walks: uno asociado al comportamiento de pasos cortos y Ć”ngulos grandes y otro asociado al comportamiento de pasos largos y los Ć”ngulos de giro pequeƱos. Como en los ejemplos anteriores, como el del willow tit, vamos a usar una variable oculta. En este caso, la variable oculta (que vamos a llamar idx), es una variable que indica si una observaciĆ³n pertenece a un comportamiento o el otro.

Ahora, en lugar de tener un Ćŗnico valor de a y b, y de rho y mu, vamos a tener uno para cada uno de los comportamientos. TambiĆ©n vamos a tener parĆ”metros (q[1, 2]) que dan la probabilidad de que el animal cambie de un comportamiento al otro. En este caso, vamos a asumir que esa probabilidad depende del comportamiento en el que se encontraba en la observaciĆ³n anterior.

cat(file = "sw.bug", "
model{
  for (t in 2:n) { 
    #Likelihood para los steps
    step[t] ~ dweib(b[idx[t]], a[idx[t]]) #b es el parƔmetro de forma
    idx[t] ~ dcat(p[t,]) # Variable oculta que identifica el comportamiento.
   
    #Likelihood para los turns (definiendo manualmente a la Wrapped-Cauchy)
    ones[t] ~ dbern( wc[t] )
    wc[t] <- (1/(2*Pi)*(1-pow(rho[idx[t]],2))/
              (1+pow(rho[idx[t]],2)-2*rho[idx[t]]*cos(turn[t]-mu[idx[t]])))/ 100
    turn[t] ~ dunif(-3.14159265359, 3.14159265359) 
    
    p[t,1] <- q[idx[t-1]] 
    p[t,2] <- 1-q[idx[t-1]] 
    }
  
  # Priors 
  a[2] ~ dunif(0, 10)
  a[1] <- a[2] + eps
  
  eps ~ dnorm(0.0, 0.01)T(0.0,) 
  
  b[1] ~ dunif(0, 8)
  b[2] ~ dunif(0, 8)
  
  mu[1] ~ dunif(-1,5)
  mu[2] ~ dunif(-1,5)
  
  rho[1] ~ dunif(0.0,1.0)
  rho[2] ~ dunif(0.0,1.0)
  
  q[1] ~ dunif(0,1)
  q[2] ~ dunif(0,1)
  
  idx[1] ~ dcat(phi[]) 
  Pi <- 3.14159265359 
}")

En el modelo, definimos la probabilidad de estar en un comportamiento en la observaciĆ³n t (p[t,1]), y como son eventos excluyentes, definimos a p[t,2] como 1 - p[t,1]. AdemĆ”s, asumimos que que la probabilidad depende del comportamiento en el que se encontraba en el tiempo anterior. Pueden ver que para el parĆ”metro de escala de la Weibull a[1] usamos un truco de sumarle una valor (eps), porque asumimos que el comportamiento \(2\) deberĆ­a tener una escala mayor que el comportamiento \(1\). Como estimamos el valor de eps a partir de los datos, si las escalas de ambos comportamientos fueran similares, el valor estimado de eps serĆ­a cercano a 0.

Antes de ajustar el modelo necesitamos definir un par de variables

idx <- numeric(n) * NA
phi = c(0.5, 0.5)  # probabilidades para el estado de la primera observaciĆ³n

ones <- numeric(n) + 1
data.sw <- list("step", "turn", "n", "phi", "ones")
inits.sw <- function() list(a = c(NA, runif(1, 0.1, 2)), b = runif(2, 1, 3), 
    mu = runif(2, 0, pi), rho = runif(2, 0, 1))
params.sw <- c("a", "b", "mu", "rho", "q", "idx")

ni <- 5000
nt <- 5
nb <- 2500
nc <- 3

sw.sim <- jags(data.sw, inits.sw, params.sw, model.file = "sw.bug", n.chains = nc, 
    n.burnin = nb, n.thin = nt, n.iter = ni)
## 
## Processing function input....... 
## 
## Done. 
##  
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 483
##    Unobserved stochastic nodes: 172
##    Total graph size: 6485
## 
## Initializing model
## 
## Adaptive phase..... 
## Adaptive phase complete 
##  
## 
##  Burn-in phase, 2500 iterations x 3 chains 
##  
## 
## Sampling from joint posterior, 2500 iterations x 3 chains 
##  
## 
## Calculating statistics....... 
## 
## Done.

Hacemos ahora el posterior predictive checking nuevamente, con los valores de las posteriores con el nuevo modelo:

n.sims <- sw.sim$mcmc.info$n.samples
ppsteps <- matrix(NA, n.sims, n)
ppturns <- matrix(NA, n.sims, n)

idx <- sw.sim$sims.list$idx
a <- sw.sim$sims.list$a
b <- sw.sim$sims.list$b
mu <- sw.sim$sims.list$mu
rho <- sw.sim$sims.list$rho

for (i in 1:n.sims) {
    ppsteps[i, ] = rweibull(n, shape = b[i, idx[i, ]], scale = (1/a[i, idx[i, 
        ]])^b[i, idx[i, ]])
    ppturns[i, ] = rwrappedcauchy(n, mu = circular(mu[i, idx[i, ]]), rho = rho[i, 
        idx[i, ]])
}

oac = acf(log(step), lag.max = (nlags - 1), plot = F)  #  ACF observada

ppac = matrix(NA, n.sims, nlags)
for (i in 1:n.sims) {
    ppac[i, ] = acf(log(ppsteps[i, ]), lag.max = (nlags - 1), plot = F)$acf
}

library(coda)
hpd <- HPDinterval(as.mcmc(ppac))

op <- par(las = 1, cex.lab = 1.5, cex.axis = 1.3)
plot(oac$lag, oac$acf, type = "b", lwd = 2, col = 2, xlab = "lag", ylab = "autocorrelation")
lines(oac$lag, hpd[, 1], lwd = 2)
lines(oac$lag, hpd[, 2], lwd = 2)

par(op)

Juan Manuel Morales. 6 de Septiembre del 2015. ƚltima actualizaciĆ³n: 2019-05-30