Uncovering ecological state dynamics with hidden Markov models

Abstract Ecological systems can often be characterised by changes among a finite set of underlying states pertaining to individuals, populations, communities or entire ecosystems through time. Owing to the inherent difficulty of empirical field studies, ecological state dynamics operating at any level of this hierarchy can often be unobservable or ‘hidden’. Ecologists must therefore often contend with incomplete or indirect observations that are somehow related to these underlying processes. By formally disentangling state and observation processes based on simple yet powerful mathematical properties that can be used to describe many ecological phenomena, hidden Markov models (HMMs) can facilitate inferences about complex system state dynamics that might otherwise be intractable. However, HMMs have only recently begun to gain traction within the broader ecological community. We provide a gentle introduction to HMMs, establish some common terminology, review the immense scope of HMMs for applied ecological research and provide a tutorial on implementation and interpretation. By illustrating how practitioners can use a simple conceptual template to customise HMMs for their specific systems of interest, revealing methodological links between existing applications, and highlighting some practical considerations and limitations of these approaches, our goal is to help establish HMMs as a fundamental inferential tool for ecologists.


Contents
Appendix A Dynamic species co-existence HMM 1 Appendix B HMM software 3 A Dynamic species co-existence HMM Here we provide additional details of an HMM formulation for species co-existence dynamics based on presence-absence data (Marescot et al., 2020). Let the states S t = A (respectively S t = B and S t = AB) indicate "site occupied by species A" (respectively by species B and by both species) and S t = U indicate "unoccupied site". Define X t,k ∈ {0, 1, 2, 3}, where 0 indicates neither species was detected, 1 indicates only species A was detected, 2 indicates only species B was detected, and 3 indicates both species were detected on the kth visit at time t. We could for example have: · · · U B AB · · · 0 · · · 0 0 · · · 2 1 · · · 3 hidden S t ∈ {AB,A,B,U} observed X t,k ∈ {0, 1, 2, 3} at multiple visits k = 1, . . . , K t − 1 t t + 1 time where ψ A (respectively ψ B ) is the probability of only species A (respectively B) being present, ψ AB is the probability of both species being present, p A (respectively p B ) is probability of detecting species A given only species A is present, r AB is the probability of detecting both species given both species are present, r Ab is the probability of detecting species A, not B, given both species are present, and r aB is the probability of detecting species B, not A, given both species are present. The state transition probability matrix Γ is composed of the following parameters: • AB is the probability that both species A and B go locally extinct between t and t + 1; • A (respectively B ) is the probability that species A goes locally extinct between t and t + 1, given both species are present at t; • ν A (respectively ν B ) is the probability that species A goes locally extinct between t and t + 1, given species B was absent at t and t + 1; • γ AB is the probability that both species A and B colonise a site between t and t + 1; • γ A (respectively γ B ) is the probability that species A colonises a site between t and t + 1, given both species are absent at t; • η A (respectively η B ) is the probability that species A colonises a site between t and t + 1, given species B was present at t and t + 1; • ω A (respectively ω B ) is the probability that species A is replaced by B between t and t + 1.

B HMM software
The computational machinery of HMMs, such as the forward and Viterbi algorithms, can be coded from scratch by a proficient statistical programmer (e.g. Zucchini et al., 2016;Louvrier et al., 2018;Santostasi et al., 2019), but recent advances in computing power and user-friendly software have made the implementation of HMMs much more feasible for practitioners. Many different HMM software packages and stand-alone programs are now available, some of which are focused on specific classes of state dynamics within the individual, population, or community level of the ecological hierarchy. However, the features and capabilities of the software are varied, and it can be challenging to determine which software may be most appropriate for a specific objective. Here we will describe some of the most popular HMM software currently available, including potential advantages and disadvantages for ecological applications.
We limit our treatment to freely available R (R Core Team, 2019) packages and stand-alone programs that we believe are most accessible to ecologists and non-statisticians.
The Comprehensive R Archive Network (https://cran.r-project.org) currently hosts 26 packages that include "hidden Markov" in their description. While most HMM packages in R include data simulation, parameter estimation, and state decoding for an arbitrary number of system states, they differ in many key respects (see Table 2  Many R packages are less general and specialise on specific HMM applications within individual-or population-level ecology. The marked package  implements many of the popular capture-recapture HMMs described in Section 3.1. Packages that spe-cialise in animal movement behaviour HMMs for telemetry data, such as those described in Section 3.1.2, include bsam (Jonsen et al., 2005), moveHMM (Michelot et al., 2016), and momentuHMM (McClintock & Michelot, 2018). The package HMMoce (Braun et al., 2018) is specifically catered for HMMs that infer location from archival tag data (e.g. light levels, depth-temperature profiles) such as those described in Section 3.1.3. Using telemetry and count data, kfdnm (Schmidt et al., 2015) can fit HMMs for population abundance and related demographic parameters such as those described in Section 3.2. The package DDD (Etienne & Haegeman, 2019) implements HMMs for macroevolutionary inference about diversification rates from phylogenetic trees such as those described in Section 3.2.2. The package openpopscr (Glennie et al., 2019) can fit spatial capture-recapture HMMs that account for unobserved animal movements when estimating population-level density and survival, such as those described in Section 3.2.3. The popular package unmarked (Fiske & Chandler, 2011) includes many of the HMMs for inferring patterns and dynamics of species occurrence from repeated presence-absence data that were described in Section 3.2.3.
There are also several stand-alone, user-friendly software programs that focus on specific HMM applications in ecology. Programs MARK (White & Burnham, 1999) and E-SURGE (Choquet et al., 2009) both provide a very general framework for implementing HMMs with individual-level capture-recapture (Pradel, 2005) or population-level presence-absence (Gimenez et al., 2014) data, including observation process error arising from non-detection (Kellner & Swihart, 2014), state uncertainty (Kendall, 2009;Kendall et al., 2012), and species misidentfication (Miller et al., 2011). Program PRESENCE (Hines, 2006)  Although not intended specifically for HMMs, it is worth noting that there are a number of 5 software programs with which these types of models can be relatively easily implemented by users with minimal statistical programming experience. For Bayesian inference using MCMC sampling (Gelman et al., 2004), these include WinBUGS/OpenBUGS (Lunn et al., 2009;Kéry & Schaub, 2011;Lunn et al., 2012), JAGS (Plummer, 2003), and Stan (Gelman et al., 2015).
There are R package interfaces for all of these programs, including R2OpenBUGS (Sturtz et al., 2005), rjags (Plummer, 2019), and rstan ( the package rbi is a complete R interface for the LibBi library (Murray, 2015). The R package TMB (Kristensen et al., 2016) generally has a steeper learning curve but can be advantageous for maximum likelihood inference (e.g. Benhaiem et al., 2018;Marescot et al., 2018), particularly for mixed HMMs that include continuous-valued random effects (Altman, 2007). From a computational point of view, neither maximum likelihood estimation nor MCMC sampling is vastly superior (cf. Patterson et al., 2017, for a more comprehensive discussion). However, MCMC samplers that include both the parameter vector (θ) and the latent states (S 1 , . . . , S T ), as commonly implemented in WinBUGS/OpenBUGS and JAGS, are inherently slow; sampling from the parameter vector only while applying the forward algorithm to evaluate the likelihood will often be preferable (Turek et al., 2016;Yackulic et al., 2020).