In TheoMichelot/hmmTMB: Fit Hidden Markov Models using Template Model Builder

knitr::opts_chunk$set(
  message = FALSE, error = FALSE, warning = FALSE,
  comment = NA
)

\tableofcontents \newpage

Introduction

In this vignette, we showcase how hmmTMB can be used to analyse time series of animal locations. This is a popular application of hidden Markov models (HMMs), to identify movement modes often interpreted as behavioural states (@patterson2009; @langrock2012). There are two specialised R packages focused on HMMs for the analysis of animal movement: moveHMM and momentuHMM (@michelot2016; @mcclintock2018). Those might be a good alternative; in particular, momentuHMM has additional functionalities for processing animal telemetry data, and to specify some complex movement models (e.g., biased random walks). The main advantage of hmmTMB over those packages is the possibility to include random effects on the model parameters, as well as non-parametric effects of covariates (with automatic smoothness selection).

The most common approach is to assume that, within each state, the animal's movement follows a correlated random walk, i.e., a model that captures the speed of movement and the degree of persistence in the movement. This is convenient because animals often alternative between periods of fast, directed movement (e.g., "transiting", "exploring"), and periods of slow, undirected movement (e.g., "resting", "foraging"). There are two main ways to model correlated random walks: based on step lengths and turning angles (@morales2004), or based on locations directly (@jonsen2005). Both approaches can be implemented in hmmTMB, and we describe them in turn below.

\tcolorbox \textbf{Data requirements}

We need to make two key assumptions about the data, which may be limiting in some animal movement applications: (a) the locations are observed at regular time intervals, and (b) the locations are observed with negligible measurement error.

In cases where these assumptions are not satisfied, another modelling approach will be needed. Various state-space models have been developed to accommodate noisy, irregular animal trajectories (@jonsen2005; @johnson2008). Relevant R packages include crawl (@johnson2008) and foieGras (@jonsen2020), and momentuHMM has dedicated functions to prepare the data using crawl before an HMM analysis. \endtcolorbox

We start by loading a few packages, including moveHMM, which contains a data set of wild haggis movement that we will use here, as well as functions for data preparation.

library(moveHMM)
library(ggplot2)
theme_set(theme_bw(13))
library(hmmTMB)

Step-and-turn model

The most common correlated random walk formulation is based on two derived variables: the step length $L_t$ (distance between successive locations) and the turning angle $\varphi_t$ (angle between successive steps). We consider the followind widely-used HMM formulation: $$ \begin{cases} L_t \mid { S_t = j } \sim \text{gamma}(\mu_j, \sigma_j)\ \varphi_t \mid { S_t = j } \sim \text{von Mises}(\theta_j, \kappa_j^2) \end{cases} $$ where, in state $j$, $\mu_j > 0$ is the mean step length, $\sigma_j$ is the standard deviation of step length, $-\pi < \theta_j \leq \pi$ is the mean turning angle, and $\kappa_j > 0$ is the concentration of turning angle.

Data preparation

In many movement studies, HMMs are used to analyse time series of step lengths and turning angles, to identify phases with different levels of speed and/or directional persistence. Those variables can directly be computed from the locations; for convenience, we use the function prepData() from moveHMM to do it here.

head(haggis_data)

# Get step lengths and turning angles
hmm_data <- prepData(haggis_data, type = "UTM")
head(hmm_data)

The columns of the data frame are:

• ID is an identifier the track or individual. This data set has three different tracks, to which we will fit a common model.

• step and angle will be the response variables for the HMM.

• x and y are the coordinates of the observed locations. They will not be used in the HMM analysis directly.

• slope and temp are environmental covariates, and we will demonstrate how their effects on the transition probabilities can be estimated.

Model specification

An HMM has two components: a hidden state process, and an observation model. Those two parts can be specified separately in hmmTMB. We first create an object of class MarkovChain for the hidden process model. Here, we just indicate that it is a 2-state process, i.e., we assume that the movement patterns are driven by two underlying states.

hid1 <- MarkovChain$new(data = hmm_data, n_states = 2)

Defining the observation process takes only a little bit more work. We need to specify the family of distribution used for each response variable; here, a gamma distribution for step length, and von Mises distribution for turning angle. Specifically, we use the gamma2 distribution, which is formulated in terms of a mean and standard deviation, rather than shape and scale.

We also need to choose initial parameter values, which will be used as a starting point for model fitting. The moveHMM package has a vignette discussing some suggested approaches to choosing those initial values: \url{https://cran.r-project.org/package=moveHMM/vignettes/moveHMM-starting-values.pdf}.

The distributions need to be stored in a list (with one element for each observed variable), and the initial parameters are in a nested list, as shown below. Finally, we can use those to create an object of class Observation.

# Observation distributions
dists <- list(step = "gamma2", angle = "vm")

# Initial observation parameters
par0 <- list(step = list(mean = c(1, 5), sd = c(1, 5)),
             angle = list(mu = c(0, 0), kappa = c(1, 5)))

obs1 <- Observation$new(data = hmm_data, dists = dists, par = par0)

Model fitting

We can now create an HMM object, which combines the hidden process and the observation model. Printing it shows the model formulation, as well as the initial parameter values that it stores. Note that, because we did not specify initial parameter values for the MarkovChain object, it defaulted to a matrix with 0.9 on the diagonal, but we could have entered it manually using the argument tpm.

hmm1 <- HMM$new(obs = obs1, hid = hid1)

hmm1

We fit the model using the fit method, which takes a few seconds. After the model is fitted, printing it shows the estimated (rather than initial) parameter values. Note that, in models where the parameters depend on covariates, only the parameters corresponding to the first row of covariates are shown, which is why the output says (t = 1). Here, the parameters are not time-varying, so this does not matter.

hmm1$fit(silent = TRUE)

hmm1

Model visualisation

The fitted model can be plotted in a few different ways. The method plot_dist creates a plot of the state-dependent distributions for the response variable specified as input. We can use such plots to help with interpretation of the states, to see how much overlap there is between distributions in different states, and to check whether the distributions seem to fit the data well.

# Plot state-dependent distributions of steps and angles
hmm1$plot_dist("step")
hmm1$plot_dist("angle")

We can also use plot_ts to show a time series plot of one or two observed variables, coloured by the Viterbi-estimated state sequence. The variables don't have to be those to which the model was fitted, so we can for example use x and y here, to create a map of the tracks. This code also highlights that the plotting functions all return ggplot objects, which can be edited as needed.

# Plot track coloured by estimated states
hmm1$plot_ts("x", "y") + 
  geom_point(size = 0.5) + 
  coord_equal()

Covariates

A popular extension of HMMs in ecology consists in including the effects of time-varying (or individual-specific) covariates on parameters of the model. This can for example be used to estimate effects of covariates on the transition probabilities of the hidden state process, and we illustrate this case here.

\tcolorbox \textbf{Aligning the covariates with the observations}

Some care is required to ensure that the covariates are correctly aligned with the response variables (e.g., step lengths and turning angles) in the input data set. In hmmTMB, the likelihood of the model is computed as $$ L_n = \bm\delta^{(1)} \bm{P}(Z_1) \bm\Gamma^{(1)} \bm{P}(Z_2) \bm\Gamma^{(2)} \cdots \bm\Gamma^{(n-1)} \bm{P}(Z_n) \bm{1}' $$ where $\bm\delta^{(1)}$ is the initial distribution of the state process, $\bm{P}(Z_t)$ is a diagonal matrix where the $j$-th element along the diagonal is $p(Z_t \vert S_t = j)$ (the likelihood of observation $Z_t$ in state $j$), $\bm\Gamma^{(t)}$ is the transition probability matrix at time $t$, and $\bm{1}'$ is a column vector of ones (@zucchini2017). That is, $\bm\Gamma^{(t)}$ is derived using the covariate values observed at time $t$.

In the movement example above, $Z_t$ is the step length between locations at times $t$ and $t+1$ (as computed by moveHMM). In this case, it makes more sense for $\bm\Gamma^{(t)}$ to be replaced by $\bm\Gamma^{(t+1)}$, then, because it represents the transition between the steps measured over the intervals $(t, t+1]$ and $(t+1, t+2]$. To achieve this, we could either shift the step lengths down by one index in the data set, or shift the covariates up by one index. Here, we follow the latter approach. \endtcolorbox

# Shift covariates up by one index to align with steps/angles
hmm_data$slope <- c(hmm_data$slope[-1], NA)
hmm_data$temp <- c(hmm_data$temp[-1], NA)

Consider that we assume that temperature has a linear effect on transition probabilities, whereas slope has an unknown non-linear effect. The latter can be specified using the syntax from the package mgcv (which implements generalised additive models in R), with the function s for a smooth term. We create a new MarkovChain, with the new formula, and the Observation object is the same as before because the observation model hasn't changed. This time, fitting the model takes a little longer because of the need to estimate the non-parametric relationship between transition probabilities and slope.

# Formula for transition probabilities (using splines)
formula <- ~ temp + s(slope, k = 8, bs = "ts")
hid2 <- MarkovChain$new(data = hmm_data, n_states = 2, formula = formula)

# Same observation model as before
obs2 <- Observation$new(data = hmm_data, dists = dists, par = par0)

# Create and fit model
hmm2 <- HMM$new(obs = obs2, hid = hid2)
hmm2$fit(silent = TRUE)

Covariate effects can be visualised with the plot method. The first argument indicates what part of the model should be plotted, and can be tpm (transition probabilities), delta (stationary state probabilities), or obspar (observation parameters). The covariates are included in the hidden process here, so we might want to visualise the first two.

hmm2$plot("tpm", var = "slope")
hmm2$plot("tpm", var = "temp")

hmm2$plot("delta", var = "slope")
hmm2$plot("delta", var = "temp")

The plots suggest there is no evidence that temperature affects transition probabilities. There seems to be a strong effect of slope, however: it is more likely to be in state 2 for intermediate values (between 10 and 30 degrees), and more likely to be in state 1 for extreme values (either smaller than 10 or larger than 30).

Position-based model

As an alternative to the model based on step lengths and turning angles, we turn to the location-based correlated random walk of @jonsen2005. If we denote as $\bm{Z}t$ the location of the animal at time $t$ (i.e., a vector of its x and y coordinate), we consider the model $$ \bm{Z}_t \mid { S_t = j } \sim N \left[ \bm{Z}{t-1} + \alpha_j (\bm{Z}{t-1}- \bm{Z}{t-2}),\ \sigma_j \bm{I}2 \right], $$ where $\bm{I}_2$ is the 2-by-2 identity matrix and, in state $j$, the parameter $\alpha_j \in \mathbb{R}$ measures the degree of autocorrelation in the location process, and $\sigma_j > 0$ measures variability. This model gives rise to correlated movement steps because, for $\alpha_j > 0$, the expected direction of change from $\bm{Z}{t-1}$ to $\bm{Z}t$ is aligned with the previous step from $\bm{Z}{t-2}$ to $\bm{Z}_{t-1}$. If $\alpha_j = 0$, this reduces to a simple Gaussian random walk where each location is normally distributed around the previous location.

Data preparation

For each coordinate, we will need to include two covariates in the observation model: the previous location ($\bm{Z}{t-1}$), and the previous step ($\bm{Z}{t-1} - \bm{Z}_{t-2}$). We derive these two variables from the data, and we remove data rows where either is missing.

for(id in unique(hmm_data$ID)) {
  ind <- which(hmm_data$ID == id)
  hmm_data$x_lag[ind] <- c(NA, hmm_data$x[ind[-length(ind)]])
  hmm_data$y_lag[ind] <- c(NA, hmm_data$y[ind[-length(ind)]])
  hmm_data$x_diff[ind] <- c(NA, NA, diff(hmm_data$x[ind])[1:(length(ind) - 2)])
  hmm_data$y_diff[ind] <- c(NA, NA, diff(hmm_data$y[ind])[1:(length(ind) - 2)])
}

hmm_data <- na.omit(hmm_data)

Model specification

Similarly to the step-and-turn analysis, we use a 2-state model. We create the hidden state model with no special arguments, but we could include covariate effects in the transition probabilities like we did above.

hid3 <- MarkovChain$new(data = hmm_data, n_states = 2)

We define the observation model as follows:

• Each observed variable (i.e., each coordinate) follows a state-specific normal distribution.

• The mean of the normal distribution depends on the previous location and on the previous step. We remove the intercept from the R formula because there is no additive constant in the mean of the normal distribution. (If there was, this would create a systematic drift towards a preferred direction.)

• We define initial values for $\sigma_1$ and $\sigma_2$, shared between the two coordinates because we assume that the movement characteristics are the same along the two dimensions. The initial mean values are not used because there is no intercept.

dists <- list(x = "norm", y = "norm")
f <- list(x = list(mean = ~ x_lag + x_diff - 1),
          y = list(mean = ~ y_lag + y_diff - 1))
par0 <- list(x = list(mean = c(0, 0), sd = c(0.3, 3)),
             y = list(mean = c(0, 0), sd = c(0.3, 3)))
obs3 <- Observation$new(data = hmm_data, dists = dists, 
                        formulas = f, par = par0)

We need to add parameter constraints to ensure that the coefficient in from of the lagged covariate $\bm{Z}{t-1}$ is not estimated (fixed to 1), and that $\alpha_j$ and $\sigma_j$ have the same value for the two observed variables (i.e., same movement parameters along x and y dimensions). We update the parameters stored in obs to set the slope parameters for $\bm{Z}{t-1}$ to 1. The function obs3$coeff_fe() can be used to remember what the order of the coefficients is. Then, we define a vector of NAs and integers that indicates which parameters are fixed (NAs), and which share a common value (integers).

obs3$update_coeff_fe(c(1, 0, 1, 0.6, log(0.3), log(3),
                      1, 0, 1, 0.6, log(0.3), log(3)))

fixpar_obs <- c(NA, 1, NA, 2, 3, 4,
                NA, 1, NA, 2, 3, 4)
names(fixpar_obs) <- rownames(obs3$coeff_fe())
fixpar <- list(obs3 = fixpar_obs)

Model fitting and results

We combine the hidden state and observation models, and pass fixpar to specify the parameter constraints, to create the HMM object.

hmm3 <- HMM$new(obs = obs3, hid = hid3, fixpar = fixpar)
hmm3$fit(silent = TRUE)

Although we can't plot the state-dependent distributions over the data like we did for step lengths and turning angles (because the distributions are covariate-dependent), we can extract the parameter estimates to interpret the states. The standard deviation parameters were $\sigma_1 = 0.74$ and $\sigma_2 = 3.18$, suggesting that there is more variability in state 2, perhaps associated with higher movement speed. The estimated persistence parameters were $\alpha_1 = -0.09$ and $\alpha_2 = 0.68$, indicating that there was virtually no persistence in state 1, and fairly strong persistence in state 2. We can also plot the trajectory, coloured by the Viterbi state sequence; from this, it is clear that state 1 captures slow clustered movement, and state 2 is fast and directed.

# Transition probability matrix
hid3$tpm()[,,1]

# Standard deviations
obs3$par()[2,,1]

# Persistence parameters
obs3$coeff_fe()[c(2,4),]

# Tracks coloured by most likely state sequence
hmm3$plot_ts("x", "y") + 
  coord_equal() + 
  geom_point(size = 0.4)

References

TheoMichelot/hmmTMB documentation built on Dec. 13, 2024, 11:52 a.m.