elk

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

Vamos a ver ejemplos de un par de modelos de random walk (RW) ajustados a datso de trayectorias registradas con GPS. El trabajo completo está aquí.

Vamos a cargar los datos y reordenarlos un poco. Separamos los "pasos" y "giros" para uno de estos elk. Los pasos son las distancias (en km) entre posiciones consecutivas de los registros de GPS (uno para cada medianoche en este caso). Nos quedamos con un vector de steps y otro de turns. Vamos a ajustar una distribución Weibull a los steps y una Wrapped Cauchy a los turns.

library(coda)
library(circular)
library(jagsUI)

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

step <- elk$km_day[2:163]
turn <- elk$turns[2:163]
n = length(step)  # número de observaciones

Como los ceros no están definidos en la Weibull, los reemplazamos por valores pequeños.

any(step == 0, na.rm = TRUE)
## [1] TRUE

step <- ifelse(step == 0, 0.01, step)

Escribimos un modelo bugs para un RW simple donde en cada movimiento elegimos 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 circular Wrapped Cauchy. El objetivo es estimar los parámetros de la Weibull y de la Wrapped Cauchy. Commo la Wrapped Cauchy no está definida en jags, tenemos que usar un truco (ones) para estimarla.

cat(file = "single.bug", 
    "
model{
  for (t in 1:n) {
    # likelihood for steps
    step[t] ~ dweib(b, a)       # b es el parámetro de forma
    # likelihood for 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)
  }
  
  # priors on movement parameters
  a ~ dnorm(0,0.001)T(0,)
  b ~ dnorm(0,0.001)T(0,)
  
  rho ~ dunif(0,1)
  mu ~ dunif(-3.1416,3.1416)
  Pi <- 3.14159265359
}"
)

Antes de ajustar el modelo tenemos que definir algunas variables

ones <- numeric(n) + 1

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

single.sim <- jags(data.single, inits.single, params.single, model.file = "single.bug", 
    n.chains = nc, n.thin = nt, n.iter = ni)
## 
## Processing function input....... 
## 
## Done. 
##  
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
##    Graph Size: 1359
## 
## Initializing model
## 
## Adaptive phase, 100 iterations x 3 chains 
## If no progress bar appears JAGS has decided not to adapt 
##  
## No burn-in specified 
##  
## Sampling from joint posterior, 2000 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,
## burn-in = 0 iterations and thin rate = 1,
## yielding 6000 total samples from the joint posterior. 
## MCMC ran for 0.13 minutes at time 2015-10-02 08:18:53.
## 
##              mean    sd     2.5%      50%    97.5% overlap0    f  Rhat
## a           1.160 0.094    0.985    1.159    1.347    FALSE 1.00 1.000
## b           0.580 0.033    0.516    0.579    0.646    FALSE 1.00 1.000
## mu          1.581 1.209   -2.672    1.852    2.968     TRUE 0.92 1.005
## rho         0.089 0.052    0.006    0.087    0.196    FALSE 1.00 1.002
## deviance 2987.042 2.677 2983.492 2986.551 2993.608    FALSE 1.00 1.000
##          n.eff
## a         3305
## b         6000
## mu        1062
## rho       1074
## deviance  6000
## 
## 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.6 and DIC = 2990.627 
## DIC is an estimate of expected predictive error (lower is better).

plot(single.sim)

Posterior Predictive Checking

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

Vamos a necesitar una función para simular valores de una Wrapped Cauchy:
rwcauchy <- function(n, mu = 0, rho = 0) {
    u = runif(n)
    V = cos(2 * pi * u)
    c = 2 * rho/(1 + rho^2)
    t <- (sign(runif(n) - 0.5) * acos((V + c)/(1 + c * V)) + mu)%%(2 * pi)
    return(t)
}

Ahora definimos cuántas simulaciones vamos a hacer y unas matrices para guardar los steps y turns simulados. Luego hacemos un loop sobre las muestras de la posterior y simulamos datos.

n.sims <- single.sim$mcmc.info$n.samples  # número de simulaciones
pps.single <- matrix(NA, n.sims, n)  # para posterior predictive steps
ppt.single <- matrix(NA, n.sims, n)  # para posterios predictive turns

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

for (i in 1:n.sims) {
    pps.single[i, ] <- rweibull(n, shape = b[i], scale = 1/a[i])  # ojo, en R usamos 1/a para la escala porque WinBUGS usa otra formulación
    ppt.single[i, ] <- rwcauchy(n, mu[i], rho[i])
}

Ahora comparamos visualmente la autocorrelación en la distancia de desplazamiento entre datos observados y datos simulados.

ppac.s = matrix(NA, n.sims, 61)
for (i in 1:n.sims) {
    ppac.s[i, ] = acf(log(pps.single[i, ]), lag.max = 60, plot = F)$acf
}

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

# matplot(oac$lag,t(ppac.s),lty=1,pch=1,col='gray')
plot(oac$lag, oac$acf, type = "b", lwd = 2, col = 2, xlab = "lag", ylab = "autocorrelation")
miquantile <- function(x) quantile(x, prob = c(0.025, 0.975))
qs = apply(ppac.s, 2, miquantile)
lines(oac$lag, qs[1, ])
lines(oac$lag, qs[2, ])

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. Vamos a ver si la cosa mejora ajustando un modelo donde el animal alterna entre dos tipos de random walks.

cat(file = "sw.bug", 
    "
model{
  for (t in 2:n) {
    # likelihood for steps
    step[t] ~ dweib(b[idx[t]], a[idx[t]]) # b is the shape parameter
    idx[t] ~ dcat(p[t,])
    # likelihood for turns
    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)
    # the probability of being in movement type 1
    p[t,1] <- q[idx[t-1]] # <- max(.000000000000,min(.999999999999, q[idx[t-1]]))
    p[t,2] <- 1-q[idx[t-1]] # <- max(.000000000000,min(.999999999999, q[idx[t-1]]))
  }
  
  # priors on movement parameters
  a[2] ~ dunif(0, 10)
  a[1] <- a[2] + eps
  
  eps ~ dnorm(0.0, 0.01)T(0.0,)
  
  b[1] ~ dunif(0.4,8)
  b[2] ~ dunif(0.4,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 
}"
)

Antes de ajustar el modelo necesitamos definir un par de variables

idx <- numeric(n) * NA  # un vector de NAs que va a contener un indicador del RW correspondiente a cada observación
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, 0.4, 3), 
    mu = runif(2, 0, pi), rho = runif(2, 0, 1))
params.sw <- c("a", "b", "mu", "rho", "q", "idx")

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

sw.sim <- jags(data.sw, inits.sw, params.sw, model.file = "sw.bug", n.chains = nc, 
    n.thin = nt, n.iter = ni)
## 
## Processing function input....... 
## 
## Done. 
##  
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
##    Graph Size: 3415
## 
## Initializing model
## 
## Adaptive phase, 100 iterations x 3 chains 
## If no progress bar appears JAGS has decided not to adapt 
##  
## No burn-in specified 
##  
## Sampling from joint posterior, 5000 iterations x 3 chains 
##  
## 
## Calculating statistics....... 
## 
## Done.

print(sw.sim)
## JAGS output for model 'sw.bug', generated by jagsUI.
## Estimates based on 3 chains of 5000 iterations,
## burn-in = 0 iterations and thin rate = 1,
## yielding 15000 total samples from the joint posterior. 
## MCMC ran for 2.588 minutes at time 2015-10-02 08:19:04.
## 
##              mean     sd     2.5%      50%    97.5% overlap0     f  Rhat
## a[1]        3.650  0.447    2.855    3.623    4.613    FALSE 1.000 1.001
## a[2]        0.186  0.109    0.045    0.161    0.456    FALSE 1.000 1.005
## b[1]        1.064  0.078    0.913    1.063    1.220    FALSE 1.000 1.000
## b[2]        1.085  0.243    0.664    1.070    1.595    FALSE 1.000 1.003
## mu[1]       2.611  0.295    2.031    2.615    3.169    FALSE 1.000 1.000
## mu[2]       0.026  0.139   -0.232    0.019    0.317     TRUE 0.561 1.000
## rho[1]      0.230  0.061    0.103    0.232    0.344    FALSE 1.000 1.000
## rho[2]      0.585  0.095    0.378    0.593    0.749    FALSE 1.000 1.000
## q[1]        0.905  0.030    0.842    0.907    0.957    FALSE 1.000 1.000
## q[2]        0.331  0.105    0.134    0.330    0.542    FALSE 1.000 1.001
## idx[1]      1.278  0.448    1.000    1.000    2.000    FALSE 1.000 1.000
## idx[2]      1.012  0.108    1.000    1.000    1.000    FALSE 1.000 1.011
## idx[3]      1.021  0.142    1.000    1.000    1.000    FALSE 1.000 1.004
## idx[4]      2.000  0.000    2.000    2.000    2.000    FALSE 1.000    NA
## idx[5]      2.000  0.000    2.000    2.000    2.000    FALSE 1.000    NA
## idx[6]      2.000  0.000    2.000    2.000    2.000    FALSE 1.000    NA
## idx[7]      2.000  0.000    2.000    2.000    2.000    FALSE 1.000    NA
## idx[8]      2.000  0.012    2.000    2.000    2.000    FALSE 1.000 1.291
## idx[9]      2.000  0.000    2.000    2.000    2.000    FALSE 1.000    NA
## idx[10]     1.229  0.420    1.000    1.000    2.000    FALSE 1.000 1.000
## idx[11]     1.038  0.192    1.000    1.000    2.000    FALSE 1.000 1.003
## idx[12]     1.034  0.182    1.000    1.000    2.000    FALSE 1.000 1.005
## idx[13]     1.999  0.023    2.000    2.000    2.000    FALSE 1.000 1.122
## idx[14]     1.953  0.212    1.000    2.000    2.000    FALSE 1.000 1.004
## idx[15]     1.078  0.269    1.000    1.000    2.000    FALSE 1.000 1.001
## idx[16]     1.016  0.125    1.000    1.000    1.000    FALSE 1.000 1.002
## idx[17]     1.007  0.085    1.000    1.000    1.000    FALSE 1.000 1.008
## idx[18]     1.034  0.180    1.000    1.000    2.000    FALSE 1.000 1.001
## idx[19]     1.012  0.108    1.000    1.000    1.000    FALSE 1.000 1.001
## idx[20]     1.001  0.035    1.000    1.000    1.000    FALSE 1.000 1.014
## idx[21]     1.002  0.042    1.000    1.000    1.000    FALSE 1.000 1.053
## idx[22]     1.001  0.036    1.000    1.000    1.000    FALSE 1.000 1.009
## idx[23]     1.004  0.059    1.000    1.000    1.000    FALSE 1.000 1.001
## idx[24]     1.023  0.150    1.000    1.000    1.000    FALSE 1.000 1.002
## idx[25]     1.895  0.306    1.000    2.000    2.000    FALSE 1.000 1.001
## idx[26]     2.000  0.012    2.000    2.000    2.000    FALSE 1.000 1.105
## idx[27]     1.203  0.402    1.000    1.000    2.000    FALSE 1.000 1.001
## idx[28]     2.000  0.000    2.000    2.000    2.000    FALSE 1.000    NA
## idx[29]     2.000  0.000    2.000    2.000    2.000    FALSE 1.000    NA
## idx[30]     2.000  0.000    2.000    2.000    2.000    FALSE 1.000    NA
## idx[31]     1.205  0.404    1.000    1.000    2.000    FALSE 1.000 1.000
## idx[32]     1.437  0.496    1.000    1.000    2.000    FALSE 1.000 1.000
## idx[33]     2.000  0.008    2.000    2.000    2.000    FALSE 1.000 1.291
## idx[34]     1.970  0.169    1.000    2.000    2.000    FALSE 1.000 1.000
## idx[35]     1.996  0.064    2.000    2.000    2.000    FALSE 1.000 1.012
## idx[36]     1.935  0.247    1.000    2.000    2.000    FALSE 1.000 1.000
## idx[37]     1.204  0.403    1.000    1.000    2.000    FALSE 1.000 1.000
## idx[38]     1.982  0.132    2.000    2.000    2.000    FALSE 1.000 1.003
## idx[39]     1.047  0.212    1.000    1.000    2.000    FALSE 1.000 1.000
## idx[40]     1.071  0.257    1.000    1.000    2.000    FALSE 1.000 1.000
## idx[41]     1.047  0.212    1.000    1.000    2.000    FALSE 1.000 1.001
## idx[42]     1.032  0.175    1.000    1.000    2.000    FALSE 1.000 1.005
## idx[43]     1.002  0.049    1.000    1.000    1.000    FALSE 1.000 1.028
## idx[44]     1.004  0.060    1.000    1.000    1.000    FALSE 1.000 1.006
## idx[45]     1.001  0.036    1.000    1.000    1.000    FALSE 1.000 1.037
## idx[46]     1.003  0.052    1.000    1.000    1.000    FALSE 1.000 1.022
## idx[47]     1.031  0.175    1.000    1.000    2.000    FALSE 1.000 1.001
## idx[48]     1.166  0.372    1.000    1.000    2.000    FALSE 1.000 1.001
## idx[49]     2.000  0.000    2.000    2.000    2.000    FALSE 1.000    NA
## idx[50]     1.677  0.468    1.000    2.000    2.000    FALSE 1.000 1.000
## idx[51]     2.000  0.000    2.000    2.000    2.000    FALSE 1.000    NA
## idx[52]     2.000  0.000    2.000    2.000    2.000    FALSE 1.000    NA
## idx[53]     2.000  0.000    2.000    2.000    2.000    FALSE 1.000    NA
## idx[54]     1.043  0.202    1.000    1.000    2.000    FALSE 1.000 1.002
## idx[55]     1.007  0.086    1.000    1.000    1.000    FALSE 1.000 1.003
## idx[56]     1.008  0.089    1.000    1.000    1.000    FALSE 1.000 1.002
## idx[57]     1.001  0.031    1.000    1.000    1.000    FALSE 1.000 1.083
## idx[58]     1.001  0.034    1.000    1.000    1.000    FALSE 1.000 1.002
## idx[59]     1.001  0.033    1.000    1.000    1.000    FALSE 1.000 1.126
## idx[60]     1.001  0.029    1.000    1.000    1.000    FALSE 1.000 1.044
## idx[61]     1.004  0.064    1.000    1.000    1.000    FALSE 1.000 1.017
## idx[62]     1.010  0.100    1.000    1.000    1.000    FALSE 1.000 1.005
## idx[63]     1.002  0.045    1.000    1.000    1.000    FALSE 1.000 1.005
## idx[64]     1.003  0.056    1.000    1.000    1.000    FALSE 1.000 1.007
## idx[65]     1.010  0.100    1.000    1.000    1.000    FALSE 1.000 1.008
## idx[66]     1.296  0.456    1.000    1.000    2.000    FALSE 1.000 1.000
## idx[67]     1.079  0.270    1.000    1.000    2.000    FALSE 1.000 1.001
## idx[68]     1.046  0.209    1.000    1.000    2.000    FALSE 1.000 1.000
## idx[69]     1.353  0.478    1.000    1.000    2.000    FALSE 1.000 1.000
## idx[70]     2.000  0.008    2.000    2.000    2.000    FALSE 1.000 1.291
## idx[71]     1.222  0.416    1.000    1.000    2.000    FALSE 1.000 1.000
## idx[72]     2.000  0.000    2.000    2.000    2.000    FALSE 1.000    NA
## idx[73]     2.000  0.000    2.000    2.000    2.000    FALSE 1.000    NA
## idx[74]     2.000  0.012    2.000    2.000    2.000    FALSE 1.000 1.291
## idx[75]     1.873  0.334    1.000    2.000    2.000    FALSE 1.000 1.001
## idx[76]     2.000  0.000    2.000    2.000    2.000    FALSE 1.000    NA
## idx[77]     1.453  0.498    1.000    1.000    2.000    FALSE 1.000 1.000
## idx[78]     2.000  0.000    2.000    2.000    2.000    FALSE 1.000    NA
## idx[79]     1.230  0.421    1.000    1.000    2.000    FALSE 1.000 1.001
## idx[80]     1.005  0.071    1.000    1.000    1.000    FALSE 1.000 1.020
## idx[81]     1.006  0.078    1.000    1.000    1.000    FALSE 1.000 1.005
## idx[82]     1.001  0.028    1.000    1.000    1.000    FALSE 1.000 1.154
## idx[83]     1.007  0.082    1.000    1.000    1.000    FALSE 1.000 1.001
## idx[84]     1.001  0.035    1.000    1.000    1.000    FALSE 1.000 1.014
## idx[85]     1.003  0.057    1.000    1.000    1.000    FALSE 1.000 1.032
## idx[86]     1.015  0.120    1.000    1.000    1.000    FALSE 1.000 1.002
## idx[87]     2.000  0.008    2.000    2.000    2.000    FALSE 1.000 1.291
## idx[88]     2.000  0.000    2.000    2.000    2.000    FALSE 1.000    NA
## idx[89]     1.015  0.122    1.000    1.000    1.000    FALSE 1.000 1.010
## idx[90]     1.003  0.053    1.000    1.000    1.000    FALSE 1.000 1.060
## idx[91]     1.002  0.039    1.000    1.000    1.000    FALSE 1.000 1.012
## idx[92]     1.053  0.225    1.000    1.000    2.000    FALSE 1.000 1.001
## idx[93]     1.117  0.321    1.000    1.000    2.000    FALSE 1.000 1.001
## idx[94]     1.004  0.060    1.000    1.000    1.000    FALSE 1.000 1.011
## idx[95]     1.002  0.050    1.000    1.000    1.000    FALSE 1.000 1.001
## idx[96]     1.001  0.031    1.000    1.000    1.000    FALSE 1.000 1.002
## idx[97]     1.003  0.053    1.000    1.000    1.000    FALSE 1.000 1.005
## idx[98]     1.004  0.062    1.000    1.000    1.000    FALSE 1.000 1.002
## idx[99]     1.002  0.044    1.000    1.000    1.000    FALSE 1.000 1.042
## idx[100]    1.001  0.032    1.000    1.000    1.000    FALSE 1.000 1.072
## idx[101]    1.002  0.041    1.000    1.000    1.000    FALSE 1.000 1.001
## idx[102]    1.003  0.054    1.000    1.000    1.000    FALSE 1.000 1.011
## idx[103]    1.001  0.031    1.000    1.000    1.000    FALSE 1.000 1.010
## idx[104]    1.001  0.024    1.000    1.000    1.000    FALSE 1.000 1.051
## idx[105]    1.001  0.036    1.000    1.000    1.000    FALSE 1.000 1.009
## idx[106]    1.001  0.033    1.000    1.000    1.000    FALSE 1.000 1.075
## idx[107]    1.008  0.090    1.000    1.000    1.000    FALSE 1.000 1.000
## idx[108]    1.005  0.074    1.000    1.000    1.000    FALSE 1.000 1.003
## idx[109]    1.001  0.036    1.000    1.000    1.000    FALSE 1.000 1.037
## idx[110]    1.001  0.028    1.000    1.000    1.000    FALSE 1.000 1.010
## idx[111]    1.029  0.168    1.000    1.000    2.000    FALSE 1.000 1.002
## idx[112]    1.019  0.137    1.000    1.000    1.000    FALSE 1.000 1.001
## idx[113]    1.004  0.059    1.000    1.000    1.000    FALSE 1.000 1.012
## idx[114]    1.024  0.153    1.000    1.000    1.000    FALSE 1.000 1.003
## idx[115]    1.012  0.107    1.000    1.000    1.000    FALSE 1.000 1.003
## idx[116]    1.007  0.081    1.000    1.000    1.000    FALSE 1.000 1.000
## idx[117]    1.001  0.035    1.000    1.000    1.000    FALSE 1.000 1.014
## idx[118]    1.001  0.029    1.000    1.000    1.000    FALSE 1.000 1.003
## idx[119]    1.007  0.083    1.000    1.000    1.000    FALSE 1.000 1.003
## idx[120]    1.001  0.036    1.000    1.000    1.000    FALSE 1.000 1.075
## idx[121]    1.002  0.044    1.000    1.000    1.000    FALSE 1.000 1.004
## idx[122]    1.006  0.076    1.000    1.000    1.000    FALSE 1.000 1.003
## idx[123]    1.005  0.073    1.000    1.000    1.000    FALSE 1.000 1.000
## idx[124]    1.001  0.038    1.000    1.000    1.000    FALSE 1.000 1.007
## idx[125]    1.004  0.064    1.000    1.000    1.000    FALSE 1.000 1.008
## idx[126]    1.018  0.133    1.000    1.000    1.000    FALSE 1.000 1.001
## idx[127]    1.006  0.076    1.000    1.000    1.000    FALSE 1.000 1.006
## idx[128]    1.001  0.035    1.000    1.000    1.000    FALSE 1.000 1.018
## idx[129]    1.021  0.143    1.000    1.000    1.000    FALSE 1.000 1.000
## idx[130]    1.004  0.063    1.000    1.000    1.000    FALSE 1.000 1.008
## idx[131]    1.001  0.027    1.000    1.000    1.000    FALSE 1.000 1.004
## idx[132]    1.015  0.122    1.000    1.000    1.000    FALSE 1.000 1.000
## idx[133]    1.003  0.055    1.000    1.000    1.000    FALSE 1.000 1.022
## idx[134]    1.007  0.085    1.000    1.000    1.000    FALSE 1.000 1.001
## idx[135]    1.004  0.061    1.000    1.000    1.000    FALSE 1.000 1.012
## idx[136]    1.003  0.053    1.000    1.000    1.000    FALSE 1.000 1.011
## idx[137]    1.002  0.045    1.000    1.000    1.000    FALSE 1.000 1.019
## idx[138]    1.003  0.051    1.000    1.000    1.000    FALSE 1.000    NA
## idx[139]    1.002  0.049    1.000    1.000    1.000    FALSE 1.000 1.015
## idx[140]    1.001  0.038    1.000    1.000    1.000    FALSE 1.000 1.013
## idx[141]    1.006  0.078    1.000    1.000    1.000    FALSE 1.000 1.001
## idx[142]    1.002  0.043    1.000    1.000    1.000    FALSE 1.000 1.012
## idx[143]    1.002  0.043    1.000    1.000    1.000    FALSE 1.000 1.019
## idx[144]    1.002  0.042    1.000    1.000    1.000    FALSE 1.000 1.050
## idx[145]    1.001  0.031    1.000    1.000    1.000    FALSE 1.000 1.002
## idx[146]    1.003  0.058    1.000    1.000    1.000    FALSE 1.000 1.027
## idx[147]    1.004  0.060    1.000    1.000    1.000    FALSE 1.000 1.007
## idx[148]    1.002  0.043    1.000    1.000    1.000    FALSE 1.000 1.012
## idx[149]    1.003  0.053    1.000    1.000    1.000    FALSE 1.000 1.006
## idx[150]    1.004  0.062    1.000    1.000    1.000    FALSE 1.000 1.001
## idx[151]    1.009  0.096    1.000    1.000    1.000    FALSE 1.000 1.004
## idx[152]    1.001  0.034    1.000    1.000    1.000    FALSE 1.000 1.002
## idx[153]    1.005  0.070    1.000    1.000    1.000    FALSE 1.000 1.002
## idx[154]    1.000  0.018    1.000    1.000    1.000    FALSE 1.000 1.071
## idx[155]    1.009  0.096    1.000    1.000    1.000    FALSE 1.000 1.003
## idx[156]    1.001  0.032    1.000    1.000    1.000    FALSE 1.000 1.007
## idx[157]    1.022  0.146    1.000    1.000    1.000    FALSE 1.000 1.001
## idx[158]    1.012  0.110    1.000    1.000    1.000    FALSE 1.000 1.000
## idx[159]    1.004  0.059    1.000    1.000    1.000    FALSE 1.000 1.001
## idx[160]    1.005  0.068    1.000    1.000    1.000    FALSE 1.000 1.004
## idx[161]    1.003  0.053    1.000    1.000    1.000    FALSE 1.000 1.023
## idx[162]    1.002  0.044    1.000    1.000    1.000    FALSE 1.000 1.016
## deviance 2755.717 16.785 2730.592 2752.897 2793.153    FALSE 1.000 1.002
##          n.eff
## a[1]      3275
## a[2]       762
## b[1]      8469
## b[2]       817
## mu[1]    15000
## mu[2]    15000
## rho[1]   14799
## rho[2]   15000
## q[1]     15000
## q[2]      3289
## idx[1]   15000
## idx[2]    3754
## idx[3]    5810
## idx[4]       1
## idx[5]       1
## idx[6]       1
## idx[7]       1
## idx[8]    7500
## idx[9]       1
## idx[10]   6881
## idx[11]   4127
## idx[12]   2898
## idx[13]   6313
## idx[14]   2532
## idx[15]   6424
## idx[16]  14400
## idx[17]   8411
## idx[18]  15000
## idx[19]  15000
## idx[20]  15000
## idx[21]   4814
## idx[22]  15000
## idx[23]  15000
## idx[24]  11744
## idx[25]   8546
## idx[26]  15000
## idx[27]   3105
## idx[28]      1
## idx[29]      1
## idx[30]      1
## idx[31]  15000
## idx[32]   6980
## idx[33]  15000
## idx[34]  15000
## idx[35]   9549
## idx[36]  15000
## idx[37]   5178
## idx[38]   8640
## idx[39]  15000
## idx[40]  15000
## idx[41]   8789
## idx[42]   3494
## idx[43]   7184
## idx[44]  15000
## idx[45]   9666
## idx[46]   8198
## idx[47]  14622
## idx[48]   4091
## idx[49]      1
## idx[50]   5519
## idx[51]      1
## idx[52]      1
## idx[53]      1
## idx[54]   5173
## idx[55]  15000
## idx[56]  15000
## idx[57]   5671
## idx[58]  15000
## idx[59]   3036
## idx[60]  12178
## idx[61]   6965
## idx[62]  10697
## idx[63]  15000
## idx[64]  15000
## idx[65]   6201
## idx[66]  15000
## idx[67]  11253
## idx[68]  15000
## idx[69]  15000
## idx[70]  15000
## idx[71]   5444
## idx[72]      1
## idx[73]      1
## idx[74]   7500
## idx[75]   6748
## idx[76]      1
## idx[77]   7955
## idx[78]      1
## idx[79]   4503
## idx[80]   4911
## idx[81]  15000
## idx[82]   3156
## idx[83]  15000
## idx[84]  15000
## idx[85]   4495
## idx[86]  15000
## idx[87]  15000
## idx[88]      1
## idx[89]   3454
## idx[90]   2670
## idx[91]  15000
## idx[92]   7742
## idx[93]   7950
## idx[94]  12269
## idx[95]  15000
## idx[96]  15000
## idx[97]  15000
## idx[98]  15000
## idx[99]   5714
## idx[100]  6245
## idx[101] 15000
## idx[102] 15000
## idx[103] 15000
## idx[104] 14992
## idx[105] 15000
## idx[106]  5576
## idx[107] 15000
## idx[108] 15000
## idx[109] 10167
## idx[110] 15000
## idx[111]  8194
## idx[112] 15000
## idx[113] 11825
## idx[114]  6716
## idx[115] 12474
## idx[116] 15000
## idx[117] 15000
## idx[118] 15000
## idx[119] 15000
## idx[120]  4473
## idx[121] 15000
## idx[122] 15000
## idx[123] 15000
## idx[124] 15000
## idx[125] 14940
## idx[126] 15000
## idx[127] 13529
## idx[128] 15000
## idx[129] 15000
## idx[130] 15000
## idx[131] 15000
## idx[132] 15000
## idx[133]  7237
## idx[134] 15000
## idx[135] 10594
## idx[136] 15000
## idx[137] 12476
## idx[138] 15000
## idx[139] 13814
## idx[140] 15000
## idx[141] 15000
## idx[142] 15000
## idx[143] 13524
## idx[144]  5334
## idx[145] 15000
## idx[146]  5295
## idx[147] 15000
## idx[148] 15000
## idx[149] 15000
## idx[150] 15000
## idx[151] 12674
## idx[152] 15000
## idx[153] 15000
## idx[154] 15000
## idx[155] 15000
## idx[156] 15000
## idx[157] 15000
## idx[158] 15000
## idx[159] 15000
## idx[160] 15000
## idx[161]  7480
## idx[162] 15000
## deviance  1558
## 
## **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 = 140.7 and DIC = 2896.428 
## DIC is an estimate of expected predictive error (lower is better).

Hacemos ahora el posterior predictive checking

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, 
        ]])  # note 1/a for scale since WinBUGS and R have different formulations
    ppturns[i, ] = rwcauchy(n, mu[i, idx[i, ]], rho[i, idx[i, ]])
}

# -- Autocorrelation in distance moved
oac = acf(log(step), lag.max = 60, plot = F)  #  ACF observada

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

# matplot(oac$lag,t(ppac),lty=1,pch=1,col='gray')
plot(oac$lag, oac$acf, type = "b", lwd = 2, col = 2)
qsw = apply(ppac, 2, miquantile)
lines(oac$lag, qsw[1, ])
lines(oac$lag, qsw[2, ])
sessionInfo()
## R version 3.2.2 (2015-08-14)
## Platform: x86_64-apple-darwin13.4.0 (64-bit)
## Running under: OS X 10.10.5 (Yosemite)
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] jagsUI_1.3.7    lattice_0.20-33 circular_0.4-7  coda_0.17-1    
## [5] knitr_1.11     
## 
## loaded via a namespace (and not attached):
##  [1] digest_0.6.8         mime_0.4             grid_3.2.2          
##  [4] knitrBootstrap_1.0.0 formatR_1.2          magrittr_1.5        
##  [7] evaluate_0.7.2       stringi_0.5-5        rjags_3-15          
## [10] boot_1.3-17          rmarkdown_0.8        tools_3.2.2         
## [13] stringr_1.0.0        markdown_0.7.7       parallel_3.2.2      
## [16] yaml_2.1.13          htmltools_0.2.6
Juan Manuel Morales. 6 de Septiembre del 2015. Última actualización: 2015-10-02