---
title: "Hidden Markov model"
author: "Georg Rüppel"
date: "`r Sys.Date()`"
output:
  rmarkdown::html_vignette:
    toc: yes
vignette: |
  %\VignetteIndexEntry{Hidden Markov model}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

Hidden Markov models (HMMs) describe the evolution of a hidden state process, $\{S_t\}$ (i.e. behavioural states), which is not directly observable through telemetry data but can be inferred from an associated sequence of observable random variables, $\{Y_t\}$, such as locations derived from telemetry.
In the case of movetrack, the user specifies the number of hidden states, $N$, which determines the range of behavioural patterns the model can represent.
These states might correspond to behaviours such as 'migratory flight', 'local movement', or 'stopover', but these are merely examples—the actual interpretation of states depends on the number chosen by the user and the data at hand.

```{r, echo=FALSE}
#| fig-cap: "Basic model structure showing the hidden state process {$S_t$} governing the real locations $z$ and the observed locations $y$ at time $t$."

DiagrammeR::grViz(" digraph {
  fontname='Helvetica,Arial,sans-serif'
  node [fontname='Helvetica,Arial,sans-serif']
  edge [fontname='Helvetica,Arial,sans-serif']

    # Weight [fillcolor='#5983b0', style=filled]

    z0 [label='z&#x209C;&#x208B;&#x2082;']
    y0 [label='y&#x2C7C;&#x208B;&#x2082;']
    z1 [label='z&#x209C;&#x208B;&#x2081;']
    y1 [label='y&#x2C7C;&#x208B;&#x2081;']
    z2 [label='z&#x209C;']
    y2 [label='y&#x2C7C;']
    z3 [label='z&#x209C;&#x208A;&#x2081;']
    y3 [label='y&#x2C7C;&#x208A;&#x2081;']
    s0 -> z0 -> y0
    s1 -> z1 -> y1
    s2 -> z2 -> y2
    s3 -> z3 -> y3

    subgraph cluster_0 {
    style=filled;
    color=lightgrey;
    node [style=filled,color=white];
        {rank=same; s0 s1 s2 s3}
    s0 [label='s&#x209C;&#x208B;&#x2082;']
        s1 [label='s&#x209C;&#x208B;&#x2081;']
        s2 [label='s&#x209C;']
        s3 [label='s&#x209C;&#x208A;&#x2081;']
        s0 -> s1 -> s2 -> s3
    label='State process';
  }
}", width = 500, height = 300)
```

In the process model, each true location $z$---that is, the actual, unobserved two-dimensional location of the animal at time $t$, governed by the hidden state process and distinct from the potentially noisy telemetry observations---is assumed to be generated by one of $N$ distributions, where $N$ is the number of states.
The hidden state process that determines which of the distributions is chosen at time $t$ is modelled as a Markov chain: the probability of being in a particular state at time $t$ depends only on the state at $t-1$.
An observation model than links these latent true locations to the coarse raw positions (see `vignette("raw_positions")`), incorporating their measurement errors obtained in the previous step.
The HMM therefore consists of three components described below.

## State process

The state process $\{S_t\}$ of a $N$-state HMM for $T$ time steps is characterised by its state transition probability matrix $\Gamma^{(t)} = (\gamma^{(t)}_{i,j})$, where $i,j = 1, \dots, N$ and $\gamma^{(t)}_{i,j} = \text{Pr}(s_{t+1} = j|s_t = i)$.
The probability of transitioning to state $s_t$ from state $s_{t-1}$ is

$$
s_t \sim \text{Categorical}(\Gamma^{(t-1)}),
    \quad 1 \leq t \leq T.
$$

## Process model

The process equation for the true locations of the animal at regular time intervals $t$, $z_t = \begin{bmatrix} z_{t, \text{lon}} \\ z_{t, \text{lat}} \end{bmatrix}$, assumes that the animal's location at time $t$ is not only dependent on the previous location, $z_{t-1}$, but also on the animal's previous movement in each coordinate, $z_{t-1} - z_{t-2}$:

$$
z_t = z_{t-1} + \lambda_n (z_{t-1} - z_{t-2}) + \epsilon_t, 
    \quad \epsilon_t \sim \text{N}(0, \Omega),
    \quad 1 \leq n \leq N,
$$

where the covariance matrix for the process variation, $\Omega$, is

$$
\Omega = 
    \begin{bmatrix} 
        \tau^2_{\epsilon, \text{lon}} & 0 \\ 
        0 & \tau^2_{\epsilon, \text{lat}} 
    \end{bmatrix}.
$$

The state-dependent parameter $\lambda_n$ ranges between 0 and 1 (i.e., $0 \leq \lambda \leq 1$) and controls how strongly the animal's current movement is influenced by its previous direction and speed.
Higher values of $\lambda_n$, approaching 1, indicate strong directional persistence---suggesting that the animal continues to move in a similar direction and at a consistent speed---typical of migratory flights.
In contrast, lower values of $\lambda_n$ reflect more random or diffusive movement patterns, where current direction and speed are less dependent on past movement---characteristic of local movements.
During stopovers, $\lambda_n$ is expected to be close to 0, reflecting that the animal did not move or remained within a certain local patch.

By default, movetrack estimates track-specific $\lambda_n$ values.
While it is also possible to use a common $\lambda_n$ across all tracks, this is not typically recommended.
A shared value for $\lambda_n$ may be advantageous in specific scenarios---such as (i) when all tracks are expected to exhibit similar movement patterns, or (ii) when individual tracks have relatively few observations, in which case pooling information across tracks can enhance model stability.
In most applications, however, allowing individual $\lambda_n$ values better captures the track-specific movement patterns of each individual.

## Observation model

The observed locations of an animal, $y_j = \begin{bmatrix} y_{j, \text{lon}} \\ y_{j, \text{lat}} \end{bmatrix}$ (i.e. the calculated coarse raw positions), often have irregular time intervals $j$, with $J$ representing the total number of observed locations.
Therefore, the true location of the animal is linearly interpolated to the time of the observation, with $w_j$ representing the proportion of the regular time interval between $t-1$ and $t$ when the observation $y_j$ was made:

$$
y_j = w_j z_t + (1 - w_j) z_{t-1} + \theta_j, 
    \quad \theta_j \sim \text{T}(0, \sigma_j), 
    \quad 1 \leq j \leq J,
$$

where $\text{T}(0, \sigma_j)$ denotes a bivariate Student's $t$-distribution with measurement error $\sigma_j = \begin{bmatrix} \sigma_{j, \text{lon}} \\ \sigma_{j, \text{lat}} \end{bmatrix}$.
The HMM is implemented using Hamilton Monte Carlo in [Stan](https://mc-stan.org).

## References

Auger-Méthé, M., Newman, K., Cole, D., Empacher, F., Gryba, R., King, A. A., ... & Thomas, L. (2021). A guide to state-space modeling of ecological time series. *Ecological Monographs*, 91(4), e01470. doi: [10.1002/ecm.1470](https://doi.org/10.1002/ecm.1470)
