Uncovering ecological state dynamics with hidden Markov models

• Ecological state inference is challenging because systems can behave nonlinearly and non-stationarily, change rapidly or subtly, and have states that are difficult or impossible to observe directly.Source: §1 · pp. 2–3
• HMMs model underlying state dynamics and observations as linked time series, disentangling them to detect state changes in noisy, time-dependent ecological phenomena.Source: §1 · p. 3
• HMM applications span identifying regime shifts, inferring population and ecosystem dynamics, detecting biodiversity from eDNA, and linking species occurrence or capture-recapture data to ecological processes.Source: §1 · p. 3
• New data streams and computing tools have fueled HMM use, but their ecological potential remains underrecognized partly because applications are framed as specialized tools rather than a general framework.Some established inferential frameworks in population, community, and movement ecology are special cases of HMMs.Source: §1 · pp. 3, 5
• The review introduces HMM structure and terminology, synthesizes applications through a flexible conceptual template, identifies promising uses and practical challenges, and provides an implementation tutorial.Its aim is to connect existing applications and establish HMMs as a fundamental tool for ecologists working with sequential data.Source: §1 · p. 5

2.1 Basic model formulation

• An HMM models sequential observations as noisy, state-dependent measurements of a finite hidden state process, with each state depending on the preceding state and observations conditionally independent given the current state.A basic model specifies initial state probabilities, transition probabilities, and state-dependent observation distributions.Source: §2.1 · pp. 7, 8, 10

2.2 Inference

• The forward algorithm computes the observation-sequence likelihood in time linear in sequence length, while maximum-likelihood or Bayesian methods estimate parameters and decoding algorithms infer states.The Viterbi algorithm finds the most likely state sequence, whereas the forward-backward algorithm estimates each time point’s state probabilities conditional on observations.Source: §2.2 · pp. 11–12

2.3 Extensions

• HMM extensions address duration dependence, ecological covariates and sequence heterogeneity, and observation processes with serial dependence or multivariate measurements.Higher-order or semi-Markov state processes can relax geometric sojourn times; covariates or random effects can represent drivers and heterogeneity, while multivariate observation models may account for dependence beyond conditional independence.Source: §2.3 · pp. 13–14

3 Ecological applications of hidden Markov models

• HMM applications are organized across ecological levels and around existential, developmental, and spatial states, with observation and hidden-state processes tailored to each system.The section uses this framework to show how HMMs can identify dynamics that are difficult or impossible to observe directly.Source: §3 · pp. 15–16

3.1.1 Existential state

• The Cormack-Jolly-Seber survival model is a two-state HMM with an absorbing dead state, estimating survival and enabling inference about when an individual most likely died.Its binary encounter data cannot distinguish death from non-detection at a given occasion; extending the model with an unborn state allows inference about birth and death using the same data.Source: §3.1.1 · pp. 17–19

3.1.2 Developmental state

• Individual developmental-state HMMs represent uncertain traits such as breeding status through state-specific survival and transition probabilities, supporting biological hypothesis tests and life-history inference.Movement HMMs infer behavioral modes from step lengths and turning angles, but their data-driven states may not correspond exactly to biologically meaningful entities.Source: §3.1.2 · pp. 20–22

3.1.3 Spatial state

• Spatial capture–recapture HMMs interpret state transitions as site-specific survival and movement, allowing decoding of likely locations during periods without observed locations.Geolocation HMMs can also infer location from indirect measurements such as light, pressure, temperature, and tidal patterns.Source: §3.1.3 · pp. 23, 25

3.2 Population level

• Population-level inference can link multiple-individual models to population concepts or model population-level observations directly without an explicit individual-process model.These approaches include linking abundance to individual survival and recruitment, or treating abundance itself as a hidden state measured by noisy counts or other observations.Source: §3.2 · p. 25, §3.2.1 · pp. 25–26

3.2.1 Existential state

• Population abundance can be derived from individual-level capture–recapture processes or represented directly as a hidden state whose transitions describe population dynamics.The direct formulation can use a Leslie matrix for transitions and count distributions such as binomial or Poisson; auxiliary data can be integrated when counts alone are insufficient.Source: §3.2.1 · pp. 25–27

3.2.2 Developmental state

• Population developmental structure—including age, sex, breeding status, fitness, and heterogeneity—can be represented by extending the alive state into a network of states and transitions.The framework treats individuals as the units driving population change, while HMM applications also extend to evolutionary processes such as those involving genes.Source: §3.2.2 · p. 27
• HMMs for molecular sequences represent evolutionary dynamics, while telemetry-based HMMs have investigated intraspecific competition, herding, and fish social behaviour (Langrock et al., 2014a).Population-level HMMs can encode age, size, sex, genotype, or metapopulation states and underlie models such as open-population N-mixture and distance-sampling models (Newman et al., 2014; Schmidt et al., 2015; Sollmann et al., 2015).Source: §3.2.2 · pp. 27–28

3.2.3 Spatial state

• Spatial capture-recapture HMMs represent population density across space and can model survival, recruitment, and movement while relating distributions to connectivity, dispersal, resources, and environmental impacts (Royle et al., 2018; Glennie et al., 2019).These models require spatial information in the data, such as spatial capture-recapture surveys.Source: §3.2.3 · p. 28

3.3 Community level

• Community-level HMM applications include systems with many species and systems with a few potentially interacting species, using multi-species presence–absence data from surveys or eDNA samples.This section introduces those two study-system types before considering existential, developmental, and spatial community states.Source: §3.3 · p. 31

3.3.1 Existential state

• Multi-species occupancy HMMs jointly infer species richness and patch occupancy, test factors shaping both and their dynamics, and can account for species never detected.These models share the HMM decomposition into initial-state, transition, and observation processes, even when fitted with hierarchical Bayesian methods.Source: §3.3.1 · p. 31

3.3.2 Developmental state

• Dynamic multi-species occupancy HMMs infer community composition, species turnover and extinction, and interspecific interactions, with co-occurrence models mostly applied to predator–prey relationships.They can also quantify α, β, and γ diversity and Hill-number measures; community assembly or succession dynamics are proposed as a potential extension, not a demonstrated application.Source: §3.3.2 · p. 32

3.3.3 Spatial state

• Dynamic multi-species occupancy and co-existence HMMs can infer changing community distributions and use local colonisation and extinction to study interspecific drivers of co-occurrence (Russell et al., 2009; Fidino et al., 2019; Marescot et al., 2019).Their more complex observation and state processes remain expressible through the HMM components δ, Γ, and P(xt).Source: §3.3.3 · pp. 32–33

3.4 Ecosystem level

• HMMs remain uncommon at ecosystem level, partly because observations are difficult to obtain and integrate across the necessary scales; a standard Markov model is an HMM assuming perfect observation.Applications include regime shifts, climate-driven dynamics, management effects, ecosystem health, precipitation, streamflow, wetlands, and land cover; more holistic integration is envisioned.Source: §3.4 · pp. 33–34

4.1 Software

• Freely available R packages vary in capabilities; depmixS4 supports multivariate sequences, covariates, constraints, and missing data, while momentuHMM also supports mixed, hierarchical, zero-inflated, and partially observed-state HMMs.The packages differ in supported features and state-dependent distributions; other R and stand-alone programs target specific ecological applications or ease implementation for users with limited statistical programming experience.Source: §4.1 · pp. 34–35, Table 2 · p. 36

4.2 Challenges and pitfalls

• HMMs suit systems whose state dynamics are Markovian and observations meet conditional-independence assumptions; residual diagnostics and ecological expertise help assess adequacy and whether richer dependence or other latent-variable models are needed.Higher-order or semi-Markov dependence can represent memory, but is more complex; other latent-variable approaches may better suit complex temporal dependence.Source: §4.2 · pp. 37–38
• Large state spaces can remain tractable with local, sparse transitions, but fully parameterized transitions may be prohibitively costly; coarsening reduces size yet can mask distinct states, so test sensitivity.Whether an HMM is feasible depends on available computing resources and state-process properties; overly coarse states can produce spurious inference.Source: §4.2 · pp. 38–39
• Model choices such as state count, covariates, and dependence are case-specific and require system expertise; some processes, like Cormack-Jolly-Seber, directly determine state count and emissions.More complex data may require experience and intuition about both the data and HMMs to identify an adequate formulation.Source: §4.2 · p. 39
• Numerical estimation can yield local maxima, label switching, or poor MCMC mixing, with risks growing as models become more complex; maximum-likelihood fitting should use many initial parameter vectors to reduce false convergence.HMM parameter estimation lacks an analytical solution, so it relies on numerical procedures.Source: §4.2 · p. 39
• Unsupervised HMM states are data-driven patterns, not necessarily biological states; locations alone cannot establish foraging, though auxiliary data or partly predefined state characteristics can improve interpretation.Labels such as “foraging” may describe an area-restricted-search pattern without proving the animal was foraging; inferred states may be crude proxies for ecological states of interest.Source: §4.2 · pp. 39–40

5 Future directions

• Ecological awareness of HMMs has increased rapidly, and their use may be even more widespread than explicit references suggest.A Web of Science search showed increasing publications and citations over the past two decades, though differing terminology and unlabelled applications may understate their prevalence.Source: §5 · p. 40, Figure 6 · p. 41
• HMMs can capture serial correlation, non-linearity, non-normality, and non-stationarity while helping ecologists adapt a conceptual template to their sequential data.The authors frame existing applications as a synthesis for linking observations to system-specific processes, rather than as isolated domain-specific examples.Source: §5 · pp. 40–41
• HMMs could integrate biotic and abiotic observations at large spatio-temporal scales to study ecosystem processes, using ecosystem-level states or lower-level dynamics within larger models.Hierarchical multi-resolution, nonparametric, and coupled HMM developments extend the capacity to represent complex data structures and relationships from individuals to ecosystems (Leos-Barajas et al., 2017a; Sherlock et al., 2013).Source: §5 · pp. 41–42
• Wider ecosystem-scale use requires familiarity with HMMs, accessible and compatible data, and discrete-state formulations with informative observation processes.Data limitations can obstruct empirical links between observation and state processes across spatio-temporal scales; ecosystem dynamics must also be faithfully conceptualized within the discrete-state framework.Source: §5 · p. 42

— Supplementary Material

• Brett T. McClintock is affiliated with NOAA’s National Marine Fisheries Service in the U.S.A.Source: — Supplementary Material · p. 68
• Roland Langrock is affiliated with Bielefeld University’s Department of Business Administration and Economics.Source: — Supplementary Material · p. 68
• Olivier Gimenez and Emmanuelle Cam are affiliated with French research institutions, including CNRS and Ifremer.Gimenez is at CNRS Centre d’Ecologie Fonctionnelle et Evolutive; Cam is at the Laboratoire des Sciences de l’Environnement Marin, Institut Universitaire Européen de la Mer, Univ. Brest, CNRS, IRD, Ifremer.Source: — Supplementary Material · p. 68
• David L. Borchers and Richard Glennie are affiliated with the University of St Andrews, while Toby A. Patterson is affiliated with CSIRO Oceans and Atmosphere in Australia.Source: — Supplementary Material · p. 68

• The hidden states are A, B, AB, and U, while each visit records whether neither, one, or both species were detected.Source: Appendix A · p. 69
• The initial state distribution assigns probabilities ψAB, ψA, and ψB to joint and single-species occupancy, with the remainder assigned to U.Source: Appendix A · pp. 69–70
• The transition matrix represents local extinction through probabilities for joint extinction, extinction given joint occupancy, and extinction conditional on the other species’ absence.Source: Appendix A · pp. 69–70
• Colonization parameters distinguish joint colonization, colonization when both species were absent, and colonization conditional on the other species’ presence; replacement probabilities are also included.Source: Appendix A · pp. 69–71
• The observation process uses species-specific detection probabilities for single-species occupancy and separate probabilities for detecting both or only one species under joint occupancy.Source: Appendix A · p. 70

Appendix B: HMM software

• R contains 26 packages listing “hidden Markov,” but many focus on categorical sequences; more flexible packages support varied state-dependent distributions and advanced HMM structures.depmixS4 supports multivariate sequences, covariates, constraints, and missing observations; momentuHMM also supports mixed and hierarchical HMMs, zero-inflated distributions, and partially observed states, but not custom-coded state-dependent distributions. Only hsmm and mhsmm are noted as supporting hidden semi-Markov models.Source: Appendix B · pp. 71–72
• Specialized R packages address ecological applications including capture-recapture, animal movement, archival-tag location, abundance, diversification, spatial capture-recapture, and species occurrence.Examples include marked, bsam, moveHMM, momentuHMM, HMMoce, kfdnm, DDD, openpopscr, and unmarked.Source: Appendix B · p. 72
• MARK and E-SURGE support capture-recapture and presence-absence HMMs with observation error, state uncertainty, and species misidentification; PRESENCE focuses on presence-absence models.RMark and RPresence provide R interfaces to MARK and PRESENCE, respectively.Source: Appendix B · p. 72
• For Bayesian or maximum-likelihood inference, users can implement HMMs with BUGS, JAGS, Stan, nimble, or TMB; jointly sampling parameters and latent states in some MCMC implementations is inherently slow.Sampling only parameters while using the forward algorithm to evaluate the likelihood is often preferable; the passage notes neither maximum likelihood nor MCMC is vastly superior overall.Source: Appendix B · p. 73

— Supplementary Tutorial

• The supplementary tutorial lists Brett T. McClintock of NOAA’s National Marine Fisheries Service among its contributors.The contributor list also names researchers affiliated with universities and research institutions in France, the United Kingdom, and Australia.Source: — Supplementary Tutorial · p. 79

S1 Introduction

• Formulating an ecologically faithful HMM, fitting it to real data, and checking model adequacy are non-trivial, so the tutorial highlights decisions and assumptions throughout the analysis.It tailors an HMM to a striated caracara’s one-hour, 1-Hz VDBA time series.Source: S1 Introduction · pp. 80–81
• The tutorial covers HMM formulation, parameter estimation, candidate-model selection, fit checking, latent-state inference, and interpretation, with an R script reproducing the caracara example in momentuHMM.The script can also be used to explore the example further.Source: S1 Introduction · p. 81

S2 Building the HMM

• An HMM models both the observation process and its underlying state process, with assumptions that contribute to its simplicity and can limit how faithfully it represents either process.The observations are realizations of the sequential process X1, . . . , XT.Source: S2 Building the HMM · p. 81

S2.1 State process

• The state sequence is a Markov chain whose next-state probabilities depend only on the current state, implying geometrically distributed state durations and a consequential memory assumption.Its evolution is determined by an N × N transition probability matrix.Source: S2.1 State process · p. 82
• Choose the number and structure of states to represent plausible ecological dynamics, balancing ecological knowledge, computational feasibility, and fit rather than relying solely on AIC or BIC (Pohle et al., 2017).In the caracara example, biologists expected four activity states: resting, minimal activity, moderate activity, and flight.Source: S2.1 State process · pp. 82–83
• Transition matrices may allow every state switch, impose biologically necessary forbidden transitions, or vary over time through covariates; initial-state probabilities may use the chain’s stationary distribution.Source: S2.1 State process · p. 83

S2.2 Observation process

• Observations are assumed conditionally independent of one another given the states, so the state process alone accounts for serial dependence.This can be unrealistic when observations remain dependent during periods spent in a single state.Source: S2.2 Observation process · p. 84
• State-dependent observation distributions can be discrete, continuous, multivariate, and covariate-dependent; positive continuous caracara VDBA was modelled with a gamma distribution.The observation model can also encode system-specific structure, such as zero detection probability in a “dead” state.Source: S2.2 Observation process · p. 84

S3.1 Overview

• HMM parameters can be estimated by direct numerical maximum likelihood, maximum likelihood via EM, or Bayesian MCMC, and efficient recursive likelihood evaluation is central to these approaches.Source: S3.1 Overview · p. 85

S3.2 Likelihood evaluation using the forward algorithm

• Under the Markov and conditional-independence assumptions, the HMM likelihood is obtained by summing over state sequences and factoring observations given their states.The state-dependent observation densities or probabilities depend on the chosen observation-process distributions.Source: S3.2 Likelihood evaluation using the… · p. 85
• The forward algorithm avoids the infeasible N^T summation over state sequences and evaluates the likelihood with cost linear in T.This can make likelihood evaluation feasible for sequences comprising millions of observations.Source: S3.2 Likelihood evaluation using the… · pp. 86–87
• The forward variable αt(j) combines the likelihood of observations through time t with information about state j at time t, and is updated step-by-step using model dependence assumptions.The update retains information about probabilities of being in the different states.Source: S3.2 Likelihood evaluation using the… · p. 86
• Starting with α1 = δP(x1), successive forward variables yield the likelihood as the sum of αT’s elements; the calculation depends on δ, Γ, and the state-dependent observation distributions.In matrix form, the likelihood multiplies the initial distribution and observation matrices with transition matrices, followed by a vector of ones.Source: S3.2 Likelihood evaluation using the… · pp. 86–87

S3.3 Model fitting

• Numerical likelihood maximisation was used to fit the VDBA HMM, with a four-state example estimating the initial distribution, transition matrix, and state-dependent distributions.Figure 2 shows estimated gamma state-dependent distributions for a striated caracara.Source: S3.3 Model fitting · p. 87, Figure 2 · p. 88
• Local likelihood maxima can be problematic in complex models, so the authors recommend trying several random initial parameter vectors during numerical optimisation.Numerical underflow is described as fairly straightforward to address.Source: S3.3 Model fitting · p. 89
• The EM algorithm is feasible and can make the maximisation problem easier than numerical likelihood maximisation in some cases, though it can be more cumbersome to implement.It can exploit recursive HMM schemes such as the forward algorithm.Source: S3.3 Model fitting · p. 89
• MCMC is a feasible fitting approach that readily enables a Bayesian HMM, using the forward algorithm to marginalise over states rather than data augmentation.The passage cites Turek et al. (2016), Patterson et al. (2017), and Yackulic et al. (2020).Source: S3.3 Model fitting · p. 89

S4 Model selection & model checking

• Candidate HMMs may vary in state number, state-dependent distribution families, or transition covariates, and their selection can follow general model-selection recommendations.The task is described as analogous to model selection in regression analysis.Source: S4 Model selection & model checking · p. 89
• Use exploratory data analysis and relevant theory to keep the candidate set small, limiting both computational demands and selection-bias risk.The passage cites Zucchini (2000) regarding selection bias.Source: S4 Model selection & model checking · p. 89
• Use model-selection criteria as guidance rather than relying on them alone.Pohle et al. (2017) is cited for a more detailed discussion of HMM model selection.Source: S4 Model selection & model checking · p. 89
• Carefully inspect promising models’ goodness of fit to determine whether they capture study-relevant data patterns and identify components needing modification.Source: S4 Model selection & model checking · p. 89

S4.1 Model selection

• Expert knowledge and the implausibility of five distinct behavioural modes in one hour motivate restricting the VDBA candidates to four or fewer states.The passage describes this as the Step 1 restriction before model selection.Source: S4.1 Model selection · pp. 89–90
• Both AIC and BIC are minimised by the VDBA model with N = 4 states.The candidate models were restricted to four or fewer states based on expert knowledge and the implausibility of five genuinely distinct behavioural modes in one hour.Source: S4.1 Model selection · pp. 89–90
• With data-driven states, AIC or BIC may support biologically unreasonable complexity because extra states can compensate for inadequate model structure, making Step 1 judgments important despite poorer fit (Leos-Barajas et al., 2017; cf. Pohle et al., 2017).Source: S4.1 Model selection · p. 90

S4.2 Model checking

• Pseudo-residuals transform each observation’s conditional distribution, given all other observations, to a standard-normal scale; deviations from normality suggest potential lack of fit.For a correctly specified model with continuous observations, the conditional probability integral transform is uniform and the resulting pseudo-residuals are standard normal.Source: S4.2 Model checking · pp. 90–91
• Pseudo-residuals should have much less correlation than the original observations, although they are not uncorrelated even under a correctly specified model.Their conditional distributions are evaluated in the context of neighbouring observations.Source: S4.2 Model checking · p. 90
• The four-state VDBA pseudo-residual QQ-plot shows tail deviations, with fewer low-acceleration observations than expected; the discrepancy is minor unless very small values matter.The QQ-plot’s dashed lines indicate point-wise 95% confidence envelopes.Source: S4.2 Model checking · p. 91, Figure 3 · p. 91
• Computing pseudo-residuals requires forward and backward variables, can be tedious, and may involve numerical underflow.HMM-specialist software can perform the calculations.Source: S4.2 Model checking · p. 91
• The fitted four-state VDBA HMM captures overall data patterns, but simulated-data checks reveal minor misfit for small observations and underestimated autocorrelation.The autocorrelation function also shows higher-than-expected correlation between neighbouring observations and underestimation at higher lags.Source: S4.2 Model checking · p. 92, Figure 4 · p. 92
• An informal alternative is to simulate observations from the fitted model and compare relevant data patterns, such as quantiles or sample autocorrelation, with the real data.Whether more complex model extensions are warranted depends on the likely importance of the misfit for the study aim.Source: S4.2 Model checking · p. 92

S5.1 Overview

• Local decoding selects the most probable state separately at each time, whereas global decoding selects the jointly most probable complete state sequence.Source: S5.1 Overview · p. 93
• Because local decoding ignores serial dependence, individually most likely states need not form the jointly most likely sequence.Source: S5.1 Overview · p. 93
• Over 95% of locally and globally decoded states agreed in the VDBA example; the choice depends on whether the question concerns a point state or the chain’s likely sequence.Figure 5 displays the globally decoded VDBA states.Source: S5.1 Overview · p. 93, Figure 5 · p. 94

S5.2 Computation

• Local decoding combines forward and backward probabilities, using αt(j)βt(j) = Pr(x1, . . . , xT, st = j), then maximizes each time point’s state posterior.The locally decoded state lt is the state value that maximizes the posterior at time t.Source: S5.2 Computation · pp. 93–94
• Global decoding uses the Viterbi recurrence qt(j) = maxi qt−1(i)γi,jf(xt | st = j), then traces backward from the maximizing state at T.The backward steps choose gt = argmaxj qt(j)γj,gt+1 for t < T.Source: S5.2 Computation · p. 94
• Only local decoding quantifies uncertainty in decoded states; Viterbi decoding returns a hard most-likely sequence without probabilities for potential misclassifications.Source: S5.2 Computation · p. 94

Interpretation

• HMM states may not match intended ecological states, so results should be interpreted cautiously and used as aids to—not substitutes for—ecological theory.Choices about state number, data streams, and HMM structure can encourage the model to reflect ecological dynamics of interest, but do not guarantee it.Source: S6 Interpretation · p. 95
• Interpretation must account for assumptions about state number, state-dependent distributions, the Markov process, and conditional independence, with model-checking deficiencies potentially limiting conclusions.Deficiencies may motivate more complex HMMs or alternative latent-variable models, or warrant restricting conclusions to robust but weak statements.Source: S6 Interpretation · p. 95
• In the VDBA example, experts associated four states with resting, minimal activity, moderate activity, and flying, assigned to 31%, 27%, 27%, and 14% of one hour’s time steps.Model checking found relatively minor issues with the marginal distribution and dependence structure; these may matter more for questions about time spent in states than broad behavioural classification.Source: S6 Interpretation · pp. 95–96
• The tutorial demonstrates an HMM workflow for uncovering ecological dynamics and highlights its computational efficiency and applicability across a wide range of ecological applications.Source: S6 Interpretation · p. 96