Movement of Elk (Cervus canadensis) reintroduced in Bancroft, Ontario
Lets look at a couple of simple Random Walk models (RW) and fit them to individual trajectories registered with GPS. More details here. First we’ll load the data for one of the tracked elk. The dataframe contains, among other things, daily steps and turning angles measured from the GPS locations. Steps are the distance in km between successive midnight locations and turns are angular differences (in radians) between successive steps. We will fit a Weibull distribution to the steps and a Wrapped Cauchy to the 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) # number of observationsAs zeros are not defined in the Weibull, we replace them with a small number.
any(step == 0, na.rm = TRUE)## [1] TRUEstep <- ifelse(step == 0, 0.01, step)Lets write a bugs model for a simple RW where for every move we choose a displacement distance from a Weibull and a turn from a Wrapped Cauchy. The goal is to estimate the parameters of the Weibull and the Wrapped Cauchy. Because the Wrapped Cauchy is not defined in JAGS we have to use a trick to simulate its pdf.
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
}"
)Before fitting the model we have to define some 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.141 minutes at time 2016-03-07 18:56:40.
##
## mean sd 2.5% 50% 97.5% overlap0 f Rhat
## a 1.160 0.097 0.977 1.156 1.353 FALSE 1.000 1.000
## b 0.579 0.033 0.517 0.579 0.645 FALSE 1.000 1.000
## mu 1.618 1.183 -2.692 1.878 2.961 TRUE 0.926 1.009
## rho 0.089 0.051 0.006 0.086 0.194 FALSE 1.000 1.001
## deviance 2987.029 2.743 2983.499 2986.499 2993.971 FALSE 1.000 1.001
## n.eff
## a 6000
## b 3436
## mu 1478
## rho 2404
## deviance 2684
##
## 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.8 and DIC = 2990.789
## DIC is an estimate of expected predictive error (lower is better).Posterior Predictive Checking
Once we fit a model it would be nice that the data we used to fit it are “possible” under the fitted model. Also, we want to see if the fitted model can reproduce some features of the data. In this case we will look at the autocorrelation between steps.
We will need a function to simulate values from a 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)
}Now we define how many simulated walks to run and set up matrices to store results. Then we write a loop to simulate random walks and store the steps and turns. To take into account the full uncertainty in out estimates, for every simulation we use values from the joint posterior distribution of the model parameters (in this case a, b, mu and rho).
n.sims <- single.sim$mcmc.info$n.samples # number of simulations
pps.single <- matrix(NA, n.sims, n) # matrix to save steps
ppt.single <- matrix(NA, n.sims, n) # matrix to save 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]) # in R we use 1/a for the scale of the Weibull because BUGS uses a different fomrulation
ppt.single[i, ] <- rwcauchy(n, mu[i], rho[i])
}Now we can take a look at the observed autocorrelation in distance moved and the one expected in the fitted model
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) # observed acf
# 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, ])
We can see that the simple RW is not cappable of reproducing the autocorrelation pattern found in the data. Maybe this pattern is due to the animal changing its movement behavior every now and then. Lets fit a model that allows that.
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
}"
)Before we fit this model we have to define some variables.
idx <- numeric(n) * NA # a vector of NAs to hold an indicator for the RW corresponding to each observation
phi = c(0.5, 0.5) # probabilities for the state of the first observation
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.627 minutes at time 2016-03-07 18:56:52.
##
## mean sd 2.5% 50% 97.5% overlap0 f Rhat
## a[1] 3.668 0.444 2.881 3.637 4.622 FALSE 1.000 1.000
## a[2] 0.191 0.114 0.047 0.164 0.476 FALSE 1.000 1.002
## b[1] 1.067 0.078 0.921 1.065 1.225 FALSE 1.000 1.000
## b[2] 1.076 0.245 0.654 1.064 1.575 FALSE 1.000 1.003
## mu[1] 2.616 0.301 2.045 2.620 3.180 FALSE 1.000 1.002
## mu[2] 0.023 0.142 -0.244 0.016 0.319 TRUE 0.551 1.001
## rho[1] 0.229 0.062 0.103 0.231 0.345 FALSE 1.000 1.000
## rho[2] 0.581 0.098 0.365 0.588 0.748 FALSE 1.000 1.001
## q[1] 0.906 0.030 0.840 0.908 0.958 FALSE 1.000 1.000
## q[2] 0.327 0.107 0.127 0.325 0.538 FALSE 1.000 1.001
## idx[1] 1.274 0.446 1.000 1.000 2.000 FALSE 1.000 1.000
## idx[2] 1.014 0.118 1.000 1.000 1.000 FALSE 1.000 1.020
## idx[3] 1.025 0.156 1.000 1.000 1.000 FALSE 1.000 1.002
## 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.016 2.000 2.000 2.000 FALSE 1.000 1.030
## idx[9] 2.000 0.000 2.000 2.000 2.000 FALSE 1.000 NA
## idx[10] 1.232 0.422 1.000 1.000 2.000 FALSE 1.000 1.001
## idx[11] 1.039 0.195 1.000 1.000 2.000 FALSE 1.000 1.002
## idx[12] 1.038 0.191 1.000 1.000 2.000 FALSE 1.000 1.000
## idx[13] 1.999 0.034 2.000 2.000 2.000 FALSE 1.000 1.079
## idx[14] 1.954 0.210 1.000 2.000 2.000 FALSE 1.000 1.000
## idx[15] 1.082 0.275 1.000 1.000 2.000 FALSE 1.000 1.000
## idx[16] 1.017 0.129 1.000 1.000 1.000 FALSE 1.000 1.005
## idx[17] 1.007 0.084 1.000 1.000 1.000 FALSE 1.000 1.000
## idx[18] 1.033 0.178 1.000 1.000 2.000 FALSE 1.000 1.004
## idx[19] 1.012 0.107 1.000 1.000 1.000 FALSE 1.000 1.002
## idx[20] 1.001 0.028 1.000 1.000 1.000 FALSE 1.000 1.039
## idx[21] 1.002 0.040 1.000 1.000 1.000 FALSE 1.000 1.018
## idx[22] 1.001 0.027 1.000 1.000 1.000 FALSE 1.000 1.099
## idx[23] 1.004 0.059 1.000 1.000 1.000 FALSE 1.000 1.006
## idx[24] 1.023 0.151 1.000 1.000 1.000 FALSE 1.000 1.001
## idx[25] 1.900 0.300 1.000 2.000 2.000 FALSE 1.000 1.000
## idx[26] 2.000 0.008 2.000 2.000 2.000 FALSE 1.000 1.291
## idx[27] 1.222 0.415 1.000 1.000 2.000 FALSE 1.000 1.000
## 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.212 0.409 1.000 1.000 2.000 FALSE 1.000 1.001
## idx[32] 1.440 0.496 1.000 1.000 2.000 FALSE 1.000 1.001
## idx[33] 2.000 0.000 2.000 2.000 2.000 FALSE 1.000 NA
## idx[34] 1.972 0.165 1.000 2.000 2.000 FALSE 1.000 1.003
## idx[35] 1.997 0.050 2.000 2.000 2.000 FALSE 1.000 1.004
## idx[36] 1.938 0.241 1.000 2.000 2.000 FALSE 1.000 1.000
## idx[37] 1.226 0.418 1.000 1.000 2.000 FALSE 1.000 1.000
## idx[38] 1.981 0.137 2.000 2.000 2.000 FALSE 1.000 1.003
## idx[39] 1.047 0.211 1.000 1.000 2.000 FALSE 1.000 1.001
## idx[40] 1.071 0.257 1.000 1.000 2.000 FALSE 1.000 1.002
## idx[41] 1.046 0.209 1.000 1.000 2.000 FALSE 1.000 1.000
## idx[42] 1.030 0.170 1.000 1.000 2.000 FALSE 1.000 1.001
## idx[43] 1.004 0.060 1.000 1.000 1.000 FALSE 1.000 1.026
## idx[44] 1.005 0.074 1.000 1.000 1.000 FALSE 1.000 1.049
## idx[45] 1.002 0.043 1.000 1.000 1.000 FALSE 1.000 1.045
## idx[46] 1.003 0.052 1.000 1.000 1.000 FALSE 1.000 1.004
## idx[47] 1.033 0.177 1.000 1.000 2.000 FALSE 1.000 1.004
## idx[48] 1.172 0.377 1.000 1.000 2.000 FALSE 1.000 1.002
## idx[49] 2.000 0.000 2.000 2.000 2.000 FALSE 1.000 NA
## idx[50] 1.680 0.466 1.000 2.000 2.000 FALSE 1.000 1.001
## 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.053 0.224 1.000 1.000 2.000 FALSE 1.000 1.000
## idx[55] 1.011 0.105 1.000 1.000 1.000 FALSE 1.000 1.007
## idx[56] 1.011 0.102 1.000 1.000 1.000 FALSE 1.000 1.007
## idx[57] 1.002 0.042 1.000 1.000 1.000 FALSE 1.000 1.088
## idx[58] 1.002 0.046 1.000 1.000 1.000 FALSE 1.000 1.119
## idx[59] 1.001 0.029 1.000 1.000 1.000 FALSE 1.000 1.081
## idx[60] 1.001 0.027 1.000 1.000 1.000 FALSE 1.000 1.004
## idx[61] 1.003 0.057 1.000 1.000 1.000 FALSE 1.000 1.015
## idx[62] 1.012 0.111 1.000 1.000 1.000 FALSE 1.000 1.001
## idx[63] 1.003 0.051 1.000 1.000 1.000 FALSE 1.000 1.012
## idx[64] 1.004 0.061 1.000 1.000 1.000 FALSE 1.000 1.002
## idx[65] 1.012 0.108 1.000 1.000 1.000 FALSE 1.000 1.001
## idx[66] 1.301 0.459 1.000 1.000 2.000 FALSE 1.000 1.000
## idx[67] 1.087 0.283 1.000 1.000 2.000 FALSE 1.000 1.000
## idx[68] 1.057 0.231 1.000 1.000 2.000 FALSE 1.000 1.000
## idx[69] 1.360 0.480 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.239 0.427 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.018 2.000 2.000 2.000 FALSE 1.000 1.189
## idx[75] 1.875 0.331 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.458 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.236 0.425 1.000 1.000 2.000 FALSE 1.000 1.001
## idx[80] 1.010 0.098 1.000 1.000 1.000 FALSE 1.000 1.008
## idx[81] 1.008 0.089 1.000 1.000 1.000 FALSE 1.000 1.001
## idx[82] 1.001 0.033 1.000 1.000 1.000 FALSE 1.000 1.075
## idx[83] 1.005 0.068 1.000 1.000 1.000 FALSE 1.000 1.017
## idx[84] 1.001 0.032 1.000 1.000 1.000 FALSE 1.000 1.019
## idx[85] 1.002 0.049 1.000 1.000 1.000 FALSE 1.000 1.015
## idx[86] 1.013 0.114 1.000 1.000 1.000 FALSE 1.000 1.001
## idx[87] 2.000 0.012 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.017 0.130 1.000 1.000 1.000 FALSE 1.000 1.001
## idx[90] 1.003 0.057 1.000 1.000 1.000 FALSE 1.000 1.006
## idx[91] 1.002 0.045 1.000 1.000 1.000 FALSE 1.000 1.007
## idx[92] 1.049 0.215 1.000 1.000 2.000 FALSE 1.000 1.001
## idx[93] 1.111 0.315 1.000 1.000 2.000 FALSE 1.000 1.000
## idx[94] 1.004 0.066 1.000 1.000 1.000 FALSE 1.000 1.015
## idx[95] 1.004 0.060 1.000 1.000 1.000 FALSE 1.000 1.001
## idx[96] 1.002 0.049 1.000 1.000 1.000 FALSE 1.000 1.005
## idx[97] 1.003 0.055 1.000 1.000 1.000 FALSE 1.000 1.002
## idx[98] 1.005 0.070 1.000 1.000 1.000 FALSE 1.000 1.000
## idx[99] 1.002 0.048 1.000 1.000 1.000 FALSE 1.000 1.012
## idx[100] 1.001 0.029 1.000 1.000 1.000 FALSE 1.000 1.044
## idx[101] 1.002 0.042 1.000 1.000 1.000 FALSE 1.000 1.026
## idx[102] 1.003 0.056 1.000 1.000 1.000 FALSE 1.000 1.024
## idx[103] 1.001 0.026 1.000 1.000 1.000 FALSE 1.000 1.059
## idx[104] 1.001 0.023 1.000 1.000 1.000 FALSE 1.000 1.088
## idx[105] 1.001 0.027 1.000 1.000 1.000 FALSE 1.000 1.060
## idx[106] 1.001 0.024 1.000 1.000 1.000 FALSE 1.000 1.051
## idx[107] 1.009 0.093 1.000 1.000 1.000 FALSE 1.000 1.003
## idx[108] 1.007 0.082 1.000 1.000 1.000 FALSE 1.000 1.002
## idx[109] 1.001 0.038 1.000 1.000 1.000 FALSE 1.000 1.013
## idx[110] 1.002 0.042 1.000 1.000 1.000 FALSE 1.000 1.018
## idx[111] 1.030 0.171 1.000 1.000 2.000 FALSE 1.000 1.003
## idx[112] 1.019 0.138 1.000 1.000 1.000 FALSE 1.000 1.010
## idx[113] 1.004 0.067 1.000 1.000 1.000 FALSE 1.000 1.010
## idx[114] 1.024 0.154 1.000 1.000 1.000 FALSE 1.000 1.001
## idx[115] 1.011 0.103 1.000 1.000 1.000 FALSE 1.000 1.001
## idx[116] 1.006 0.077 1.000 1.000 1.000 FALSE 1.000 1.013
## idx[117] 1.001 0.036 1.000 1.000 1.000 FALSE 1.000 1.069
## idx[118] 1.001 0.034 1.000 1.000 1.000 FALSE 1.000 1.099
## idx[119] 1.006 0.078 1.000 1.000 1.000 FALSE 1.000 1.007
## idx[120] 1.001 0.036 1.000 1.000 1.000 FALSE 1.000 1.017
## idx[121] 1.002 0.039 1.000 1.000 1.000 FALSE 1.000 1.006
## idx[122] 1.006 0.077 1.000 1.000 1.000 FALSE 1.000 1.000
## idx[123] 1.006 0.080 1.000 1.000 1.000 FALSE 1.000 1.003
## idx[124] 1.002 0.043 1.000 1.000 1.000 FALSE 1.000 1.037
## idx[125] 1.004 0.061 1.000 1.000 1.000 FALSE 1.000 1.023
## idx[126] 1.017 0.131 1.000 1.000 1.000 FALSE 1.000 1.002
## idx[127] 1.009 0.093 1.000 1.000 1.000 FALSE 1.000 1.005
## idx[128] 1.001 0.037 1.000 1.000 1.000 FALSE 1.000 1.038
## idx[129] 1.019 0.138 1.000 1.000 1.000 FALSE 1.000 1.003
## idx[130] 1.006 0.075 1.000 1.000 1.000 FALSE 1.000 1.004
## idx[131] 1.001 0.032 1.000 1.000 1.000 FALSE 1.000 1.007
## idx[132] 1.013 0.114 1.000 1.000 1.000 FALSE 1.000 1.001
## idx[133] 1.004 0.064 1.000 1.000 1.000 FALSE 1.000 1.021
## idx[134] 1.007 0.086 1.000 1.000 1.000 FALSE 1.000 1.007
## idx[135] 1.004 0.062 1.000 1.000 1.000 FALSE 1.000 1.015
## idx[136] 1.004 0.063 1.000 1.000 1.000 FALSE 1.000 1.025
## idx[137] 1.002 0.042 1.000 1.000 1.000 FALSE 1.000 1.005
## idx[138] 1.002 0.047 1.000 1.000 1.000 FALSE 1.000 1.028
## idx[139] 1.003 0.050 1.000 1.000 1.000 FALSE 1.000 1.002
## idx[140] 1.001 0.037 1.000 1.000 1.000 FALSE 1.000 1.010
## idx[141] 1.006 0.075 1.000 1.000 1.000 FALSE 1.000 1.002
## idx[142] 1.002 0.039 1.000 1.000 1.000 FALSE 1.000 1.028
## idx[143] 1.002 0.039 1.000 1.000 1.000 FALSE 1.000 1.012
## idx[144] 1.002 0.041 1.000 1.000 1.000 FALSE 1.000 1.005
## idx[145] 1.001 0.036 1.000 1.000 1.000 FALSE 1.000 1.037
## idx[146] 1.003 0.059 1.000 1.000 1.000 FALSE 1.000 1.014
## idx[147] 1.004 0.060 1.000 1.000 1.000 FALSE 1.000 1.007
## idx[148] 1.002 0.039 1.000 1.000 1.000 FALSE 1.000 1.025
## idx[149] 1.002 0.045 1.000 1.000 1.000 FALSE 1.000 1.002
## idx[150] 1.003 0.056 1.000 1.000 1.000 FALSE 1.000 1.015
## idx[151] 1.007 0.082 1.000 1.000 1.000 FALSE 1.000 1.002
## idx[152] 1.001 0.031 1.000 1.000 1.000 FALSE 1.000 1.002
## idx[153] 1.004 0.061 1.000 1.000 1.000 FALSE 1.000 1.000
## idx[154] 1.001 0.031 1.000 1.000 1.000 FALSE 1.000 1.010
## idx[155] 1.009 0.095 1.000 1.000 1.000 FALSE 1.000 1.011
## idx[156] 1.001 0.028 1.000 1.000 1.000 FALSE 1.000 1.106
## idx[157] 1.022 0.147 1.000 1.000 1.000 FALSE 1.000 1.007
## idx[158] 1.012 0.109 1.000 1.000 1.000 FALSE 1.000 1.006
## idx[159] 1.002 0.045 1.000 1.000 1.000 FALSE 1.000 1.002
## idx[160] 1.004 0.064 1.000 1.000 1.000 FALSE 1.000 1.004
## idx[161] 1.004 0.064 1.000 1.000 1.000 FALSE 1.000 1.009
## idx[162] 1.002 0.045 1.000 1.000 1.000 FALSE 1.000 1.022
## deviance 2756.511 17.854 2730.399 2753.507 2797.429 FALSE 1.000 1.001
## n.eff
## a[1] 14298
## a[2] 1090
## b[1] 15000
## b[2] 690
## mu[1] 5611
## mu[2] 2753
## rho[1] 4020
## rho[2] 3537
## q[1] 8054
## q[2] 1678
## idx[1] 15000
## idx[2] 1728
## idx[3] 8711
## idx[4] 1
## idx[5] 1
## idx[6] 1
## idx[7] 1
## idx[8] 15000
## idx[9] 1
## idx[10] 3302
## idx[11] 5497
## idx[12] 15000
## idx[13] 4899
## idx[14] 15000
## idx[15] 15000
## idx[16] 6360
## idx[17] 15000
## idx[18] 3798
## idx[19] 15000
## idx[20] 14989
## idx[21] 15000
## idx[22] 5890
## idx[23] 15000
## idx[24] 15000
## idx[25] 15000
## idx[26] 15000
## idx[27] 13906
## idx[28] 1
## idx[29] 1
## idx[30] 1
## idx[31] 3692
## idx[32] 1933
## idx[33] 1
## idx[34] 6403
## idx[35] 15000
## idx[36] 15000
## idx[37] 6919
## idx[38] 8485
## idx[39] 9445
## idx[40] 3635
## idx[41] 15000
## idx[42] 15000
## idx[43] 5175
## idx[44] 1749
## idx[45] 5517
## idx[46] 15000
## idx[47] 3620
## idx[48] 2115
## idx[49] 1
## idx[50] 1886
## idx[51] 1
## idx[52] 1
## idx[53] 1
## idx[54] 15000
## idx[55] 6348
## idx[56] 6344
## idx[57] 2751
## idx[58] 1641
## idx[59] 6286
## idx[60] 15000
## idx[61] 9640
## idx[62] 15000
## idx[63] 15000
## idx[64] 15000
## idx[65] 15000
## idx[66] 15000
## idx[67] 13596
## idx[68] 15000
## idx[69] 8772
## idx[70] 15000
## idx[71] 4873
## idx[72] 1
## idx[73] 1
## idx[74] 5768
## idx[75] 4470
## idx[76] 1
## idx[77] 4009
## idx[78] 1
## idx[79] 4482
## idx[80] 6436
## idx[81] 15000
## idx[82] 5576
## idx[83] 6077
## idx[84] 15000
## idx[85] 13814
## idx[86] 15000
## idx[87] 7500
## idx[88] 1
## idx[89] 15000
## idx[90] 15000
## idx[91] 15000
## idx[92] 10517
## idx[93] 10577
## idx[94] 7300
## idx[95] 15000
## idx[96] 15000
## idx[97] 15000
## idx[98] 15000
## idx[99] 15000
## idx[100] 12178
## idx[101] 10367
## idx[102] 6449
## idx[103] 11532
## idx[104] 9227
## idx[105] 10306
## idx[106] 14992
## idx[107] 15000
## idx[108] 15000
## idx[109] 15000
## idx[110] 15000
## idx[111] 4742
## idx[112] 2485
## idx[113] 10996
## idx[114] 15000
## idx[115] 15000
## idx[116] 6117
## idx[117] 4912
## idx[118] 3803
## idx[119] 10748
## idx[120] 15000
## idx[121] 15000
## idx[122] 15000
## idx[123] 15000
## idx[124] 6873
## idx[125] 5543
## idx[126] 15000
## idx[127] 10681
## idx[128] 8739
## idx[129] 7632
## idx[130] 15000
## idx[131] 15000
## idx[132] 15000
## idx[133] 5805
## idx[134] 9213
## idx[135] 8667
## idx[136] 4871
## idx[137] 15000
## idx[138] 7841
## idx[139] 15000
## idx[140] 15000
## idx[141] 15000
## idx[142] 11113
## idx[143] 15000
## idx[144] 15000
## idx[145] 10167
## idx[146] 10229
## idx[147] 15000
## idx[148] 12304
## idx[149] 15000
## idx[150] 10490
## idx[151] 15000
## idx[152] 15000
## idx[153] 15000
## idx[154] 15000
## idx[155] 5054
## idx[156] 4997
## idx[157] 3128
## idx[158] 7431
## idx[159] 15000
## idx[160] 15000
## idx[161] 13602
## idx[162] 10793
## deviance 1515
##
## **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 = 159.2 and DIC = 2915.714
## DIC is an estimate of expected predictive error (lower is better).Lets perform a posterior predictive check for this model
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) # observed acf
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, ])
Juan Manuel Morales. 6 de Septiembre del 2015. Última actualización: 2016-03-07