I’m slowly building up a model on based on GPS collar data for a small forest carnivore. The data consist of a time series of GPS fixes and accelerometer readings collected every fifteen minutes or half hour. All GPS fixes have some amount of measurement error, but because the animal is in the forest some fixes have more error than others, up to the limit of having no fix if the GPS was not in view of any satellites. We do however, always have an accelerometer measurement. One of the goals of our research is to identify when the animal is resting and the location of resting sites (often structures like tree cavities), but also to generally characterize the animals movement. Because these animals rest in cover, we tend to have more failed GPS fixes or high measurement error fixes when we presume the animal is resting in a resting site. Eventually, I want to build up to hidden Markov model with measurement error. However, for now I’m trying to understand the measurement error piece with a simple correlated random walk with measurement error.
Using some Stan and Rcode that a colleague obtained from a recent workshop on animal movement at ISEC (which I believe Juan Morales wrote) as a basis, I have the following R and Stan code for simulating and fitting a correlated random walk with measurement error. The issue is that for this relatively simple model, I am having trouble getting the model to converge without divergent transitions, exceeding the treedepth, or getting a warning about BFMI. I was eventually able to eliminate divergent transitions and treedepth exceedances by setting adapt_delta to 0.975 and max_treedepth to 20. But with 4000 iterations over 4 chains and a 10 hour runtime I still get BFMI warnings. This is for a simple model with only 200 data points in the time-series. That doesn’t give me a lot of hope for scaling this approach up to account for totally missing GPS fixes and state switching.
The model is able to return the provided parameter values even with divergent transition warnings. Indeed the fit never looks bad (except for the warnings). But I’m wondering if there is some inherent pathology here that I’m missing? Is there a better way to parameterize this?
Here is the R code to simulate:
library(CircStats)
library(rstan)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
T <- 200
scale <- 1
shape <- 1.5
m <- 0
rho <- 0.8
locs <- matrix(0, T, 2)
steps <- rep(NA, T - 1)
turn_angle <- rep(NA, T - 1)
turn_angle[1] <- runif(1,-pi, pi)
for(t in 2:T){
steps[t - 1] <- rweibull(1, shape = shape, scale = scale)
if(t > 2)
turn_angle[t - 1] <- rwrpcauchy(1, location = turn_angle[t - 2] + m + pi, rho = rho) - pi
locs[t, 1] <- locs[t - 1, 1] + cos(turn_angle[t - 1]) * steps[t - 1]
locs[t, 2] <- locs[t - 1, 2] + sin(turn_angle[t - 1]) * steps[t - 1]
}
plot(locs[,1],locs[,2], type = "l", asp = 1, xlab="", ylab="", lwd=1)
s_obs <- 0.8
xo <- locs[,1] + rnorm(dim(locs)[1], 0, s_obs)
yo <- locs[,2] + rnorm(dim(locs)[1], 0, s_obs)
plot(xo,yo, type = "l", asp = 1, xlab="", ylab="", lwd=1)
data_list <- list(
T, xo, yo
)
params <- c("m", "rho", "scale", "shape", "x", "y",
"sigma", "steps", "turn_angle")
fit <- stan(file="crw_oe.stan", data=data_list, iter = 4000,
chains = 4, pars = params)
And Stan code to fit:
functions{
//Function that returns the log pdf of the wrapped-Cauchy
real wrappedCauchy_lpdf(real y, real rho, real mu) {
return(log(1/(2*pi())) + log((1-rho^2)/(1+rho^2-2*rho*cos(y-mu))));
}
}
data {
int<lower=0> T; // length of the time series
vector[T] xo;
vector[T] yo;
}
parameters {
real<lower=0, upper=1> rho;
real<lower=-pi(), upper=pi()> m;
real<lower=0> shape;
real<lower=0> scale;
real<lower=0> sigma;
vector<lower=0>[T - 1] steps;
vector[T - 1] turn_angle;
}
transformed parameters{
vector[T] x;
vector[T] y;
x[1] = 0;
y[1] = 0;
for(t in 2:T){
x[t] = x[t - 1] + cos(turn_angle[t - 1]) * steps[t - 1];
y[t] = y[t - 1] + sin(turn_angle[t - 1]) * steps[t - 1];
}
}
model {
turn_angle[1] ~ uniform(-pi(), pi());
m ~ normal(0,0.5);
rho ~ beta(4,1);
scale ~ normal(0,1);
shape ~ gamma(2,1);
sigma ~ normal(0, .25);
for(t in 2:T){
steps[t - 1] ~ weibull(shape, scale);
if (t > 2)
turn_angle[t - 1] ~ wrappedCauchy(rho, turn_angle[t - 2] + m);
}
for (t in 1:T) {
xo[t] ~ normal(x[t], sigma);
yo[t] ~ normal(y[t], sigma);
}
}
