Hierarchical continuous time trend model for animal movement

[DEHewitt]

HI @nicholasjclark,

Picking up from our conversation in #39 - I am wanting to account for temporal autocorrelation in a time series of crab movements. We expect that these crabs are moving downstream to spawn in the ocean and want to model the influence of water temperature and the lunar cycle on this.

Each individual was tagged with an acoustic transmitter within a river system and monitored using fixed-position receivers. The tags transmit randomly every 2-5 minutes and the receivers have an approximate detection radius of 200-300 m (i.e., 0.1 km² listening area), and for the most part these crabs aren't fast swimmers. Based on the timeseries of detections I have classified movements as (where y in data denotes the movements):

1. Stationary: detection at t is on the same receiver as the detection at t - 1 (y == 0)
2. Downstream: detection at t is on a receiver that is downstream of the detection at t - 1 (y == 1)
3. Upstream: detection at t is on a receiver that is upstream of the detection at t - 1 (very rare, n = 66 out of ~13k observations; y == 0)

We tagged crabs in two different size classes (i.e., below or above their legal harvest size) and over 4 different tagging trips (Nov, Jan, Feb, Mar). Where I see an issue with autocorrelation arising is that once an individual started moving downstream they did so relatively rapidly (for a crab). So, there is an initial y == 1 when the crab switched to migratory behaviour, which we hypothesize is driven by covariates at that time, followed by a series of other y == 1 that may still depend on those covariates but also depend on the fact that the crab has 'decided' to migrate. If I were to fit a model with mgcv, it'd look like this:

m <- gam(
  formula = y ~
    # Size class intercept to set up by-factor smooths
    size_class + 
    
    # Water temperature
    s(temp, bs = "tp", by = size_class) +
    
    # Lunar phase (in radians) by size class
    s(lunar, bs = "cc", by = size_class) + 
    
    # When they were tagged
    #s(cohort, bs = "re") +
    
    # Individual-level variation nested within each cohort
    s(transmitter, bs = "re", by = cohort),
  data = data,
  family = binomial(link = "logit"),
  method = "REML",
  select = TRUE,
  knots = list(lunar = c(0, 2 * pi))
)

The other (secondary) issue is that since the crabs move slow relative to the detection area of the receivers is that sometimes a crab will be moving downstream (y == 1), but it will record multiple detections on the same receiver (y == 0) since it doesn't get out of the detection range quick enough - this strikes me as an observation model issue... maybe it could also be dealt with by doing some data-wrangling / thinning (but I am trying to adopt an ethos of finding a model to fit the data rather than making the data fit the model). In either case, one issue at a time may be wisest.

Thanks for taking the time to read this and help (if possible).

Cheers,
Dan

[nicholasjclark]

Hi @DEHewitt, thanks for that helpful description. My main question is:

If an animal switches to a movement state (y = 1), does it tend to stay 1 for the remaining duration? In other words, are they switching from a 'non-migratory' to a 'migratory' state based on some external cues and then they migrate to the ocean?

If that is the case, it sounds like perhaps you are more interested in some kind of Hidden Markov Model to try and estimate when each crab has 'switched' to a migratory state. Would this be correct, or am I on the wrong track here?

[DEHewitt]

Hi @nicholasjclark,

I think you're definitely on the right track. A hidden Markov model would potentially be good here, but I don't have the typical metrics from these tracks to fit one (e.g., step length, turning angle). The array of receivers is somewhat dispersed across a constricted river system, so even calling them tracks is generous. The data is essentially a capture-recapture timeseries (but with continuous recapture effort). I think the most appropriate analysis is probably a multistate capture-recapture model, but at this stage my coauthors prefer that I pursue a GLM/GAM. Is the kind of thing I described originally possible with your package? FWIW I am hoping to convince them that we go down the capture-recapture path, but am interested in if mvgam can be useful for modelling animal tracking timeseries in this way.

[nicholasjclark]

Thanks @DEHewitt. This is a really interesting problem. Looking at your second point:

The other (secondary) issue is that since the crabs move slow relative to the detection area of the receivers is that sometimes a crab will be moving downstream (y == 1), but it will record multiple detections on the same receiver (y == 0) since it doesn't get out of the detection range quick enough - this strikes me as an observation model issue

Makes me think that you could perhaps view this from the lens of an occupancy model, where the latent binary state is 'migrating' vs 'not migrating' and you have some detection error. Given that the crabs must all move downstream, detection error could be driven by the distance among receivers (you are more likely to mis-detect a migrating crab if there is a bigger distance from the current receiver to the next one, for example). The environmental covariates and other dynamic processes could then be used in the process to model inform predictions of whether or not the crab was migrating at a given time. I believe you could do this in {flocker}, which would still allow you to use GAMs for the relevant model subcomponents but maybe get more explicit recognition of the detection error process. Currently, models with binary latent states are not available in {mvgam}, though they will be added in the near future. But if you don't want to do this, then the best approach in my package may be to let each crab have its own CAR(1) process and use the hierarchical GAM components as you have mentioned above to partially pool their environmental responses. This is a bit more computationally demanding, especially if you have many crabs, but it might be the only way to adequately account for the uneven spacing of detections. I could help you to modify the model file so that the CAR(1) $\rho$ and $\sigma$ parameters can be learned hierarchically, which would make sense in this case I believe.