WAIC or LOO-CV for a Hidden Markov Model of animal movement

[Ben Lambert]

Hi,

I have estimated a Hidden Markov Model (HMM) for animal movement states (with the help of Bob! See the attached.) The marginalisation of the latent states is handled via dynamic programming using the forwards/backwards algorithm.

I have a number of model permutations which I would like to compare using some measure of predictive performance. For example, I would like to compare a two-state HMM with a three-state one.

My question is - does anyone have any ideas how to estimate WAIC (or LOO-CV) without repeatedly partitioning the data into T-1 points and a single datum? The model takes about half an hour to estimate for ~18,000 data points, and so doing this for all points is clearly not feasible. 

Best,

Ben

[Jonah Gabry]

You could maybe start with one of the variants of K-fold cross validation, which would definitely be more computationally appealing.

Have you checked what happens if you do the PSIS-LOO approximation in the loo package? I'm not sure that's what you want for this model, but I'd be curious to know what happens.

Jonah

[Bob Carpenter]

I think with a single time series, you need to
fit the model on observations y[1] to y[t-1]
and look at the predicted log density for
y[t].

With multiple time series (as in multiple animals),
you might be able to cross-validate across them.

I thought you and Aki and Andrew would know what
to do.

- Bob

> --
> You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

[Jonah Sol Gabry]

On Fri, Jul 8, 2016 at 12:57 AM, Bob Carpenter <ca...@alias-i.com> wrote:

I think with a single time series, you need to
fit the model on observations y[1] to y[t-1]
and look at the predicted log density for
y[t].

Well, LOO can be appropriate for time series, but not always. Depends on what you're focused on, e.g. it's ok if you are considering the ability of the model to describe the structure in the observed time series. But if you're focused on prediction for new data, then yes some form of sequential cross-validation seems ideal (but I thought too computationally demanding in this case. How big is T?). If you do try sequential cross-validation then I would probably condition on at least several data points, i.e. don't start with t = 2.  

[Ben Lambert]

Hello all,

Thanks for your messages. Just trying to understand - what do you mean by sequential cross-validation?

I was thinking that I could use a variant of k-folds for time series where I essentially do the following:

1 | 2

1 2 | 3

1 2 3 | 4

1 2 3 4 | 5

where the | indicates a boundary between a training and testing set. However, to me it is a bit unclear how to do this with a HMM. The code below shows how I estimate the model. Basically - to estimate the states - we find the most probable state for each time period, determined by considering the previous time periods' best state.

So apart from the actually fold method used, my question is how to do prediction, and how to estimate the log probability of this?

Best,

Ben

functions{

  ## Function that returns the log pdf of the wrapped-Cauchy

  real wrappedCauchy(real aPhi, real aRho, real aMu) {

    return(- log(2*pi()) + log((1-aRho^2)/(1+aRho^2-2*aRho*cos(aPhi-aMu))));

  }

  

  ## Function that returns a random sample from wrapped-Cauchy by first

  ## sampling from the Cauchy, then using modular arithmetic

  real wrappedCauchy_rng(real aRho, real aMu){

    real aCauchyRand;

    real aWCRand;

    

    aCauchyRand = cauchy_rng(pi()+aMu,-log(aRho));

    if (aCauchyRand>0){

      aWCRand = fmod(aCauchyRand,2*pi()) - pi();

    } else{

      aWCRand = -fmod(-aCauchyRand,2*pi()) + pi();   

    }

    return(aWCRand);

  }

  

}

data {

  int<lower=0> N;

  vector<lower=-pi(), upper=pi()>[N] phi;

  vector<lower=0>[N] l;

  int<lower=1> K;

}

parameters {

  simplex[K] theta[K];

  vector<lower=-pi(), upper=pi()>[K] mu;

  vector<lower=0,upper=1>[K] rho;

  positive_ordered[K] alpha;

}  

model {

  vector[K] log_theta[K];

  vector[K] lp;

  for (k in 1:K) log_theta[k] = log(theta[k]);

  lp = rep_vector(-log(K), K); ## Initial state log-probability

  for (n in 1:N) {

    for (k in 1:K)

      lp[k]

        = log_sum_exp(log_theta[k] + lp)

        + wrappedCauchy(phi[n], rho[k], mu[k])

        + exponential_lpdf(l[n] | alpha[k]);

  }

  target += log_sum_exp(lp);

  alpha[1] ~ normal(0.001,0.0001);

  alpha[2] ~ normal(0.005,0.0001);

  

}

generated quantities {

  int<lower=1,upper=K> y_star[N];

  real log_p_y_star;

  {

    int back_ptr[N, K];

    real best_logp[N, K];

    real best_total_logp;

    for (k in 1:K)

      best_logp[1, K] = wrappedCauchy(phi[1], rho[k], mu[k])

                      + exponential_lpdf(l[1] | alpha[k]);

    for (t in 2:N) {

      for (k in 1:K) {

      best_logp[t, k] = negative_infinity();

        for (j in 1:K) {

          real logp;

          logp = best_logp[t-1, j]

          + log(theta[j, k]) + wrappedCauchy(phi[t], rho[k], mu[k])

                      + exponential_lpdf(l[t] | alpha[k]);

          if (logp > best_logp[t, k]) {

            back_ptr[t, k] = j;

            best_logp[t, k] = logp;

          }

        }

      }

    }

    log_p_y_star = max(best_logp[N]);

    for (k in 1:K)

      if (best_logp[N, k] == log_p_y_star)

      y_star[N] = k;

      for (t in 1:(N - 1))

      y_star[N - t] = back_ptr[N - t + 1,

      y_star[N - t + 1]];

  }

}

[Krzysztof Sakrejda]

On Thursday, July 7, 2016 at 3:24:19 PM UTC-4, Ben Lambert wrote:

Hi,

I have estimated a Hidden Markov Model (HMM) for animal movement states (with the help of Bob! See the attached.) The marginalisation of the latent states is handled via dynamic programming using the forwards/backwards algorithm.

I have a number of model permutations which I would like to compare using some measure of predictive performance. For example, I would like to compare a two-state HMM with a three-state one.

Just a word of caution about predictions in

this context.  Lots of time-series models will produce decent

single time-step predictions just because they are set up as

stepwise models.  To really see if you have a decent model 

you need multi-step predictions on a meaningful time scale.

In other words, I can tell you where that elk is going to be

fifteen minutes after the last GPS fix, it'll be near where you

last saw it, and if you can't beat that model (y[t]=y[t-1]) an 

hour or two out you might want to rethink your modeling approach.

With infectious disease count time series we've been taking data

like this and partitioning it into meaningful blocks to do cross validation (since

dropping a single observation doesn't make much sense anyway)

but we've also taken an iterative approach where for (almost) each

t* you predict y[t>t*] conditional on y[t<=t*].  This mimics how data

would arrive and lets you look at prediction quality as a function

of t-t*.  It's a good comparison to make against models like AR(1) 

and y[t]=y[t-1] to see if you can do better.

Krzysztof

[Jonah Sol Gabry]

Re: sequential cross-validation, both what you described and what Krzysztof suggested (predict y[t>t*] conditional on y[t<=t*]) would fall under that umbrella. I think Krzysztof has much more experience with the kind of model you're fitting than I do, so I would trust his advice here over mine. 

Jonah

P.S. After talking about these HMMs for animal movement with Bob recently I do plan to spend some time on these at some point. It would be nice to have an implementation in rstanarm, as I guess these models are getting more and more popular. 

--
You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/Nb_yXztKVFY/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+...@googlegroups.com.

[Bob Carpenter]

> On Jul 8, 2016, at 10:28 AM, Krzysztof Sakrejda <krzysztof...@gmail.com> wrote:

...

> Just a word of caution about predictions in
> this context. Lots of time-series models will produce decent
> single time-step predictions just because they are set up as
> stepwise models. To really see if you have a decent model
> you need multi-step predictions on a meaningful time scale.

Do you mean predicting the entire rest of the series, like

p(y[(t+1):T] | y[1:t])

?

> In other words, I can tell you where that elk is going to be
> fifteen minutes after the last GPS fix, it'll be near where you
> last saw it, and if you can't beat that model (y[t]=y[t-1]) an
> hour or two out you might want to rethink your modeling approach.

It's a little more subtle than that. The elk tend to either
mill about, or walk in a straightish line. So it's not just
a random diffusion around where they're at for each time step.

But then when you zoom out scales, it does look like a diffusion
again. Without longer-scale history data, say for fixed migrations
paths, or covariates, such as location of food and water and terrain
features, it's hard to tell where an animal's going to go.

The real problem is that y[t] = y[t-1] isn't a model. It doesn't provide
a density for p(y[t] | y[t-1]). So how do I evaluate against a non-model?
With mean square error or something? In the mixture HMM case, there are two
components, change in heading and distance moved---together with the current
bearing and state (transiting, foraging, sleeping, etc.), they give you a
prediction of where you'll be next.

- Bob

[Krzysztof Sakrejda]

On Friday, July 8, 2016 at 1:02:23 PM UTC-4, Bob Carpenter wrote:

> On Jul 8, 2016, at 10:28 AM, Krzysztof Sakrejda <krzysztof...@gmail.com> wrote:

...

> Just a word of caution about predictions in
> this context.  Lots of time-series models will produce decent
> single time-step predictions just because they are set up as
> stepwise models.  To really see if you have a decent model
> you need multi-step predictions on a meaningful time scale.

Do you mean predicting the entire rest of the series, like

  p(y[(t+1):T] | y[1:t])

?

Yes.

 

> In other words, I can tell you where that elk is going to be
> fifteen minutes after the last GPS fix, it'll be near where you
> last saw it, and if you can't beat that model (y[t]=y[t-1]) an
> hour or two out you might want to rethink your modeling approach.

It's a little more subtle than that.  The elk tend to either
mill about, or walk in a straightish line.  So it's not just
a random diffusion around where they're at for each time step.

Yes, so after a basic exploratory look at the data you are suggesting 

that a two-state model of this sort will do better than random diffusion

at predicting animal movement.  Fine, I'm saying test that. Additionally, 

if you're going to talk about modeling animal movement test it more

than one time step ahead and see how far ahead your generative model

gets you. In general you'll do better than diffusion on short term

time scales and then worse again and the real question is how big that

window is.

  

But then when you zoom out scales, it does look like a diffusion
again.   Without longer-scale history data, say for fixed migrations
paths, or covariates, such as location of food and water and terrain
features, it's hard to tell where an animal's going to go.

Exactly, this is a really cool area that has a lot of potential for combining this

GPS/satellite data with more detailed behavioral work.  The difficulty

with it is that it's being driven by tech (GPS/satellite tags) and

animal movement modeling is at diminishing returns with that.  The missing piece

is the less flashy behavioral work (as you alude to).  

 

The real problem is that y[t] = y[t-1] isn't a model.  It doesn't provide
a density for p(y[t] | y[t-1]).  So how do I evaluate against a non-model?
With mean square error or something?  In the mixture HMM case, there are two
components, change in heading and distance moved---together with the current
bearing and state (transiting, foraging, sleeping, etc.), they give you a
prediction of where you'll be next.

Sure it's a point prediction.  You can use any error you want, for infectious disease

we've done mean absolute error and relative mean absolute error  (MAE1/MAEref)

where MAEref is for some reference model (sometimes as simple as y[t]=y[t-1]). Typically

we evaluate these as a function of prediction horizon or something like that.

Krzysztof

 

- Bob

[Bob Carpenter]

> On Jul 8, 2016, at 2:17 PM, Krzysztof Sakrejda <krzysztof...@gmail.com> wrote:
>
>
>
> On Friday, July 8, 2016 at 1:02:23 PM UTC-4, Bob Carpenter wrote:
>
> > On Jul 8, 2016, at 10:28 AM, Krzysztof Sakrejda <krzysztof...@gmail.com> wrote:
>
> ...
>
> > Just a word of caution about predictions in
> > this context. Lots of time-series models will produce decent
> > single time-step predictions just because they are set up as
> > stepwise models. To really see if you have a decent model
> > you need multi-step predictions on a meaningful time scale.
>
> Do you mean predicting the entire rest of the series, like
>
> p(y[(t+1):T] | y[1:t])
>
> ?
>
> Yes.
>
>
> > In other words, I can tell you where that elk is going to be
> > fifteen minutes after the last GPS fix, it'll be near where you
> > last saw it, and if you can't beat that model (y[t]=y[t-1]) an
> > hour or two out you might want to rethink your modeling approach.
>
> It's a little more subtle than that. The elk tend to either
> mill about, or walk in a straightish line. So it's not just
> a random diffusion around where they're at for each time step.
>
>
> Yes, so after a basic exploratory look at the data you are suggesting
> that a two-state model of this sort will do better than random diffusion
> at predicting animal movement. Fine, I'm saying test that.

That is, of course, what the ecologists do. In lots of ways,
but mainly with BIC.

So what I'm asking is, what's the baseline? It can't just be
a point mass. I suppose it could be a uniform turning angle
and then some fitted distance moved?

> Additionally,
> if you're going to talk about modeling animal movement test it more
> than one time step ahead and see how far ahead your generative model
> gets you. In general you'll do better than diffusion on short term
> time scales and then worse again and the real question is how big that
> window is.

But I'd think if the model's reasonable, it'll do better at longer
time scales, too, than a simpler model. Even though it won't be able to
make an accurate prediction, it should be able to estimate the density
of everywhere it might be. The correlation among bearings will mean
that you can explain longer walks ---it's a kind of overdispersion.

> The real problem is that y[t] = y[t-1] isn't a model. It doesn't provide
> a density for p(y[t] | y[t-1]). So how do I evaluate against a non-model?
> With mean square error or something? In the mixture HMM case, there are two
> components, change in heading and distance moved---together with the current
> bearing and state (transiting, foraging, sleeping, etc.), they give you a
> prediction of where you'll be next.
>
> Sure it's a point prediction. You can use any error you want,

Well, I want to use log loss to compare probabilistic predictions.
That is, I want to lay out an area and ask what the probability
is that the animal will be in it. Or at least evaluate the density
of the point the animal's actually at.

> for infectious disease
> we've done mean absolute error and relative mean absolute error (MAE1/MAEref)
> where MAEref is for some reference model (sometimes as simple as y[t]=y[t-1]). Typically
> we evaluate these as a function of prediction horizon or something like that.

But I want to do full Bayes for my predictions. I'd integrate
the error measure over the posterior predictive to derive an
expected error and bounds:

INT err(y[(t+1):T], y'[(t+1):T]) p(y'[(t+1):T] | y[1:t]) dy'.

- Bob

[Krzysztof Sakrejda]

If your more complex model does better than a simpler model at

shorter and longer time scales I think that's a desirable property but

you don't get it for free.  In ecology there's a nice series of articles about

this between Sugihara and Hartig that's unfortunately couched as

parametric versus non-parametric models.  

 

Even though it won't be able to
make an accurate prediction, it should be able to estimate the density
of everywhere it might be.  The correlation among bearings will mean
that you can explain longer walks ---it's a kind of overdispersion.  

> The real problem is that y[t] = y[t-1] isn't a model.  It doesn't provide
> a density for p(y[t] | y[t-1]).  So how do I evaluate against a non-model?
> With mean square error or something?  In the mixture HMM case, there are two
> components, change in heading and distance moved---together with the current
> bearing and state (transiting, foraging, sleeping, etc.), they give you a
> prediction of where you'll be next.
>
> Sure it's a point prediction.  You can use any error you want,

Well, I want to use log loss to compare probabilistic predictions.
That is, I want to lay out an area and ask what the probability
is that the animal will be in it.  Or at least evaluate the density
of the point the animal's actually at.

Yeah, fair enough, for that particular approach something like AR(1)

(for infectious disease) or diffusion (for animal movement) would 

be a better comparison.  For longer time scales I don't think these

HMM's are going to do much better than diffusion.  The length of

the window where they do better is what's interesting.

 

> for infectious disease
> we've done mean absolute error and relative mean absolute error  (MAE1/MAEref)
> where MAEref is for some reference model (sometimes as simple as y[t]=y[t-1]). Typically
> we evaluate these as a function of prediction horizon or something like that.

But I want to do full Bayes for my predictions.  I'd integrate
the error measure over the posterior predictive to derive an
expected error and bounds:

INT err(y[(t+1):T], y'[(t+1):T]) p(y'[(t+1):T] | y[1:t]) dy'.

I think we're on the same page here, except that you're really against doing

comparisons to point mass models because you can't do this integration.

That's fine, there are other time series reference models that you can do 

this integration with (AR(1) or diffusion).  My main point was that it's almost always

 interesting to do the comparison and to also check how far into the future

a fancy model does better than a reference model.  

Krzysztof

 

- Bob

[Bob Carpenter]

Thanks for the pointer (and for focusing on evaluation).
I'll have to check those out. This topic never came up at ISEC,
but it should have been central to all these papers!

> Well, I want to use log loss to compare probabilistic predictions.
> That is, I want to lay out an area and ask what the probability
> is that the animal will be in it. Or at least evaluate the density
> of the point the animal's actually at.
>
> Yeah, fair enough, for that particular approach something like AR(1)
> (for infectious disease) or diffusion (for animal movement) would
> be a better comparison. For longer time scales I don't think these
> HMM's are going to do much better than diffusion. The length of
> the window where they do better is what's interesting.

I think you're right, because on a long enough time scale, a
mixture of resting, foraging, and directed movement is going to
look like a diffusion of some sort. At shorter scales, though,
they can manage some level of directedness to movement and account
for the mixture of consistent shorter-range and longer-range
movements. Most of the focus on longer-range movement seemed to
be on covariates (where is there food, where are there hills, where
other fish or whales or panthers are) and constraints (dolphins
can't move across land), or they just directly estimated something
like a point process.

> > for infectious disease
> > we've done mean absolute error and relative mean absolute error (MAE1/MAEref)
> > where MAEref is for some reference model (sometimes as simple as y[t]=y[t-1]). Typically
> > we evaluate these as a function of prediction horizon or something like that.
>
> But I want to do full Bayes for my predictions. I'd integrate
> the error measure over the posterior predictive to derive an
> expected error and bounds:
>
> INT err(y[(t+1):T], y'[(t+1):T]) p(y'[(t+1):T] | y[1:t]) dy'.
>
> I think we're on the same page here, except that you're really against doing
> comparisons to point mass models because you can't do this integration.

I'm OK with it in other cases, I just don't know how to do
it in these diffusive continous processes. It's not like
we're estimating a simple parameter and can use mean square error
or something. At least I don't know how to.

> That's fine, there are other time series reference models that you can do
> this integration with (AR(1) or diffusion). My main point was that it's almost always
> interesting to do the comparison and to also check how far into the future
> a fancy model does better than a reference model.

Very good point. They should've been doing this in all those
ISEC papers. I'm still not 100% sure how to do this because I don't
have longer term data and I don't know much about diffusion processes.
But I'd love to see an example so I could learn how to do it.

- Bob

[Aki Vehtari]

On Friday, July 8, 2016 at 2:19:09 PM UTC+3, Ben Lambert wrote:

I was thinking that I could use a variant of k-folds for time series where I essentially do the following:

1 | 2

1 2 | 3

1 2 3 | 4

1 2 3 4 | 5

Yes, you could do this and if these correspond to blocks of data then you will be able to analyse the multistep prediction performance up to the size of the block size.

Krzysztof proposes 

1 | 2 3 4 5

1 2 | 3 4 5

1 2 3 | 4 5

1 2 3 4 | 5

 

which may be useful if the model is able to predict longer than one block.

You could also use "ghosting" by Marshall and Spiegelhalter to improve stability if removing just one block does not affect much the hyperparameter posteriori.

where the | indicates a boundary between a training and testing set. However, to me it is a bit unclear how to do this with a HMM. The code below shows how I estimate the model. Basically - to estimate the states - we find the most probable state for each time period, determined by considering the previous time periods' best state.

 Since you select just the most probable state, you can't use use importance sampling LOO or WAIC and thus you have to re-run the inference and selection with different data divisions. Blocking will help to reduce the amount of different data divisions,

Aki