Selecting the Number of States in Hidden Markov Models - Pitfalls, Practical Challenges and Pragmatic Solutions

1.8 MB·

by  Jennifer Pohle

Selecting the Number of States in Hidden Markov Models — Pitfalls,
Practical Challenges and Pragmatic Solutions
Jennifer Pohle1∗, Roland Langrock1, Floris M. van Beest2∗ and Niels Martin Schmidt2∗
7
1 1Bielefeld University, Germany
0
2 2Aarhus University, Denmark

n
a
J
0 Abstract
3
We discuss the notorious problem of order selection in hidden Markov models, i.e. of selecting an

]
E adequatenumberofstates,highlightingtypicalpitfallsandpracticalchallengesarisingwhenanalyzingreal
M data. Extensive simulations are used to demonstrate the reasons that render order selection particularly
. challenginginpracticedespitetheconceptualsimplicityofthetask. Inparticular,wedemonstratewhywell-
t
a
t established formal procedures for model selection, such as those based on standard information criteria,
s
[ tend to favor models with undesirably large numbers of states. We also offer a pragmatic step-by-step

1 approach together with comprehensive advice for how practitioners can implement order selection. Our
v
3 proposed strategy is illustrated with a real-data case study on muskox movement.
7
6
Keywords: animal movement, information criteria, Markov chain, selection bias, unsupervised learning
8
0
.
1
0 1 INTRODUCTION
7
1
: Hidden Markov models (HMMs) are flexible time series models for sequences of observations that are driven
v
i
X by underlying, serially correlated state sequences. Originating from the domain of speech recognition, they
r have found successful applications in various areas such as robotics, finance, economics and social science (cf.
a
Chapter 14 of Zucchini et al., 2016). Over the last couple of years, they have also emerged as an increasingly
popular statistical tool for the analysis of ecological time series data, where they have proven to be natural
statistical models for animal movement data (Patterson et al., 2016), general animal behavior data (DeRuiter
etal.,2016),capture-recapturedata(Pradel,2005),anddistancesamplingsubjecttoavailabilitybias(Borchers
et al., 2013), to name but a few.
In this paper, we discuss order selection in (finite-state) HMMs, i.e. how to select the number of states
of an HMM. While conceptually this appears to be a simple model selection task, in practice it remains a
notoriously difficult challenge faced by practitioners. Whether or not order selection involves difficulties very
∗Correspondingauthor,[email protected]
1
much depends on the purpose of an HMM-based analysis. We distinguish three main types of applications
of HMMs: forecasting, classification (in a supervised learning context), and general inference on the data-
generating process (unsupervised learning). Order selection is most challenging in the latter case, and we
therefore focus on this application specifically in this paper.
In principle, when a maximum likelihood approach is taken it is conceptually straightforward to use in-
formation criteria, e.g. the Akaike Information Criterion (AIC), the Bayesian Information Criterion (BIC), or
variations thereof, to select between models with different numbers of states (cf. Celeux and Durand, 2008;
Bacciet al.,2014,Zucchiniet al.,2016). However,especiallyintheecologicalliteratureithasbeenclaimed—
and sometimes demonstrated — that traditional model selection criteria, and especially the AIC, often lead
totheselectionofmuchlargernumbersofstatesthanexpecteda priori(Deanet al.,2012;vandeKerket al.,
2015; Langrock et al., 2015; DeRuiter et al., 2016; Li and Bolker, in press).
Thepreferencetoincludemanystates,particularlyinecologicalsettings,cantosomeextentbeexplainedby
the complexity of such data sets. In addition to the features that actually motivate the use of state-switching
models, such as multimodality and autocorrelation (which together are often indicative of the presence of
correlatedlatentstates),realdataoftenexhibitfurtherstructure,suchasoutliers,seasonalityorheterogeneity
between individuals. When neglecting these features in the formulation of an HMM, then additional states
maybeabletocapturethisignoreddatastructureandthereforeprovideabettermodelfitthanmodelswitha
lower, but more realistic number of hidden states (Langrock et al., 2015; DeRuiter et al., 2016; Li and Bolker,
in press). For analyses where the interest lies on the interpretation of the states, or on the general dynamics
of the state process, this behavior of traditional model selection criteria (e.g. AIC, BIC) is highly undesirable.
As an example, consider the case of animal movement data. In this context, the states can intuitively be
interpreted as proxies for the behavioral states of an animal (e.g. resting, foraging or traveling), and primary
interest usually focuses on identifying the internal and external drivers of behavioral processes. Therefore,
in the ideal case, an HMM applied to an animal’s movement data can yield a deeper understanding of the
behaviorofsaidanimal. However, asoutlinedabove, traditionalmodelselectioncriteriaoftenpointtomodels
with large numbers of states which, crucially, may not be biologically interpretable anymore.
Itisourviewthat,todate,thereisnowidelyacceptedprocedurethatcouldresolvethepracticalproblems
described above related to order selection in HMMs. In this paper, we claim that it is probably not possible
to design a one-size-fits-all formal selection method that would reliably yield sensible numbers of states in
practice. Therefore, in this paper, we suggest a more pragmatic approach to choose the number of states in
practical applications, which takes into account formal criteria for guidance, but also stresses the importance
of the study aim, expert knowledge and model checking procedures.
The paper is organized as follows. In Section 2, we lay out the HMM basics, including a brief review of
the types of ways in which HMMs are applied, and provide an overview of theoretical and practical aspects
related to order selection. In Section 3, we use simulation studies to demonstrate how additional states in an
HMM can capture neglected structure inthe data, leading standard information criteria to often overestimate
the true number of hidden states. In Section 4, we discuss how to pragmatically choose an adequate number
2
of states, and provide practical advice and guidance. Section 5 gives a real-data case study with muskox
movement data, illustrating how to implement our pragmatic approach to order selection in HMMs.
2 HIDDEN MARKOV MODEL BASICS
2.1 BASIC FORMULATION OF HMMS
An HMM is a doubly stochastic process in discrete time, the structure of which, in its most basic form, is
displayedinFigure1. Themodelassumestheobservations, X t=1,2,...,T ,tobedrivenbyanunderlying
t
{ | }
(unobserved) state sequence, S t=1,2,...,T . In movement ecology, the states are usually interpreted as
t
{ | }
proxies of the behavioral states of the animal observed (Patterson et al., 2016). The state process is usually
assumed to be an N-state Markov chain, such that
Pr(S S ,S ,...,S )=Pr(S S ).
t t−1 t−2 1 t t−1
| |
Thus, given the present state, future states are independent of past states. Without loss of generality,
throughout the paper we additionally assume the Markov chain to be stationary, unless explicitly stated
(cid:0) (cid:1)
otherwise. The state-switching probabilities are summarized in the transition probability matrix Γ = γ ,
ij
where γ =Pr(S =j S =i), i,j =1,...,N.
ij t t−1
|
Xt 1 0Xt 0 Xt+1 (observed)
−
... St 1 0St 0 St+1 ... (hidden)
−
Figure 1: Dependence structure of an HMM in its most basic form.
In addition to the Markov property, it is usually assumed that the observations are conditionally indepen-
dent of each other, and of past states, given the current state:
p(X X ,X ,..,X ,S ,S ,...,S )=p(X S ).
t t−1 t−2 1 t t−1 1 t t
| |
Here p is used as a general symbol to denote either a probability mass function (if X is discrete-valued) or a
t
densityfunction(ifX iscontinuous-valued). Thus,thedistributionofeachobservedvariableX ,t=1,...,T,
t t
is completely determined by the current state S . Together with the Markov property, this assumption gives
t
the basic model structure as directed acyclic graph depicted in Figure 1.
3
In practice, especially with animal behavior data, the observed process will often be multivariate, i.e.
X = (X ,...,X ). In that case, a commonly made additional assumption is that the m variables are
t 1t mt
also conditionally independent, given the current state. This is the so-called contemporaneous conditional
independence assumption:
m
(cid:89)
p(X S )= p(X S ).
t t it t
| |
i=1
In this basic form, the HMM formulation is completed by choosing the number of hidden states, N, and the
class(es)ofstate-dependentdistributions(alsocalledemissiondistributions). Inpractice,themodelparameters
will then need to be estimated, which is usually accomplished using numerical maximum likelihood (Zucchini
et al., 2016). In the next section, we discuss in which ways HMMs can be used in applied statistics, and any
implications of the particular type of application on the issue of choosing the number of states, N.
2.2 TYPES OF HMM-BASED ANALYSES
Asindicatedintheintroduction,therearethreemaintypesofwaysinwhichHMMscanbeapplied: forecasting,
classification and general inference. (One could argue that a possible fourth type of application is as a tool to
make likelihood inference tractable, as demonstrated for example in Pedersen et al., 2011, or in Langrock and
King, 2013 — but that type of application is not relevant for the present work.)
First, HMMs can be used for forecasting future values of the observed time series, as typically done
in econometric time series analysis (see, e.g., Hamilton and Susmel, 1994). In these instances, HMMs are
not usually regarded as good representations of the true data-generating process, and instead are typically
used merely as tools to accommodate features of the observed time series that may be difficult to capture
otherwise, using standard time series models (e.g. ARIMA- or GARCH-type models). The main purpose of
this application of HMMs is usually not to learn something about the data-generating process, but instead to
accurately predict future observations. A specific example is given by the application of HMMs to financial
share returns (Bulla, 2011): here it is not actually reasonable to assume that there is a finite number of
(market) states, yet the corresponding models can produce good forecast distributions.
Second, in the machine learning literature, HMMs are often used for (state) classification in a supervised
learning context (e.g. for speech recognition, gesture recognition, activity recognition, etc.). In these settings,
an HMM is trained using data where the underlying states are known, and subsequently applied to new,
unlabeled data with the aim to recognize the underlying states. In ecology, HMMs are sometimes applied in
this way in order to detect animal behavior states (Broekhuis et al., 2014; see also the discussion in Leos-
Barajas et al., 2016). However, these applications are relatively rare in ecology, since training data, where the
states are directly observed, need to be available, which is usually difficult to realize in the field. In this type
of application of HMMs, the choice of the number of states is not an issue, as the states and their meaning
are predefined.
Third, in an unsupervised context, HMMs are used to learn something about the data-generating process,
without defining the role of the states a priori. Especially in movement ecology, this is the standard way in
4
which HMMs are applied, with the aim of inferring novel aspects related to the behavioral process (Morales
et al., 2004; Patterson et al., 2009; Langrock et al., 2012). While practitioners may have some expectations
regarding the number of states also in the unsupervised context, the identification of the true, or at least a
suitable number of states in general still remains a primary aim of empirical studies. Thus, the unsupervised
learning context is where order selection in HMMs constitutes the biggest challenge, so it is this case that we
focus on in the current paper.
2.3 MODEL SELECTION FOR HMMS — THEORY AND PRACTICE
In practical applications of HMMs, users need to at least (i) specify the dependence assumptions made within
themodel(typicallytheMarkovpropertyandconditionalindependenceoftheobservations, giventhestates),
(ii) decide on the class of distributions used for the state-dependent process (e.g. normal distributions), and
(iii)selectthenumberofstates, N. Inaddition, itmaybenecessaryto(iv)decidewhichcovariatestoinclude
the model.
It is our experience that in most practical applications of HMMs, model selection focuses on (iii) and, if
applicable, (iv), with (i) and (ii) specified with little or no investigation into the corresponding goodness-of-
fit of the resulting models (though there are of course exceptions to this). For the model selection involved
in both (iii) and (iv), when a maximum likelihood approach is taken, then information criteria such as the
AIC ( 2log +2p, where is the maximum value of the likelihood function and p is the number of model
− L L
parameters) or the BIC ( 2log +plog(T)) are typically used.
− L
Similarly as the BIC, the integrated completed likelihood (ICL) criterion proposed by Biernacki et al.
(2001) takes into account model evidence, but additionally considers the relevance of partitions of the data
into distinct states, as obtained under the model. The ICL criterion approximates the integrated complete-
data likelihood, i.e. the joint likelihood of the observed values x=(x ,...,x ) and its associated underlying
1 T
state sequence s = (s ,...,s ) using a BIC-like approximation. As the true state sequence is unknown, it
1 T
is replaced by the Viterbi-decoded state sequence sˆ, i.e. the most probable state sequence under the model
considered. With (x,sˆ) denoting the (approximate) complete-data likelihood, the ICL criterion is defined
c
L
as
ICL= 2log (x,sˆ)+plog(T).
c
− L
As in case of the AIC and the BIC, the model is chosen that leads to the smallest value of the criterion.
In the context of HMMs, the simulation studies provided by Celeux and Durand (2008) indicate that ICL
may actually underestimate the number of states of the HMM in certain scenarios. This can be explained
by the preference of the ICL criterion for models with non-overlapping classes (i.e. states in the context of
HMMs). Despite its intuitive appeal in HMM-based clustering tasks, the ICL has not yet widely been used
by practitioners working with HMMs, such that experience regarding its practical use is limited.
Cross-validated likelihood, as suggested in Celeux and Durand (2008), constitutes another alternative,
data-driven approach to model selection, which focuses mostly on predictive performance. In the simulation
5
studiesprovidedinCeleuxandDurand(2008),cross-validationdoesnotseemtooutperformBICinitsability
tofindanadequatenumberofstates. Additionally,cross-validationmethodscanbecomeverycomputationally
intensive, which becomes particularly problematic with the increasingly large telemetry data sets collected in
movement ecology, where the estimation of a single model can easily take several hours.
From the theoretical perspective, the behavior of any of these criteria, and of ICL in particular, is still
poorly understood. The purpose of this paper is not to contribute to the theoretical development of order
selection methods in HMMs, but rather to discuss practical challenges faced by practitioners and how to
address them. For an overview of the state of the art of the theoretical development of model selection in
HMMs, see Chapter 15 in Capp´e et al. (2005). Additional formal criteria, largely similar in spirit to those
detailed above, are discussed in Celeux and Durand (2008) and Bacci et al. (2014).
In practice, to select the number of states, HMMs with different numbers of states are typically fitted to
the data until the AIC/BIC/ICL has reached its lowest value. A key assumption of both AIC and BIC is that
the true data-generating process can indeed be represented by one of the candidate models considered. If the
candidate models are at least good approximations of the true data-generating process, then this procedure
can still be expected to work well (Zucchini, 2000). However, when working with complex ecological data,
neitheranHMMnorinfactanyothercomputationallyfeasiblestatisticalmodelcanbeexpectedtobeagood
approximation of the true process, which will typically involve various complex patterns. Thus, intuitively it
is clear that undesirable behavior of say AIC or BIC may occur in practical applications of HMMs to complex
patterns, andindeedthisiswhathasbeenobservedinvariousempiricalapplicationsinthepast, inparticular
regardingorderselection(see,e.g.,Deanet al.,2012,vandeKerket al.,2015,Langrocket al.,2015,DeRuiter
et al., 2016). The purpose of this paper is to shed some light on why order selection using these standard
criteria is doomed to fail in highly complex empirical applications, and to discuss pragmatic solutions to deal
with this challenge.
2.4 RUNNING EXAMPLE — ANIMAL MOVEMENT
Toillustrateboththeproblemsoforderselectionandtheirpragmaticsolution,weusetheillustratingexample
of animal movement modeling throughout this article. In general, it is our view that an intuitive running
example helps to highlight the practical issues arising in order selection. We focus on the particular example
of animal movement as it constitutes one of the most natural and intuitive applications of HMMs, yet at the
same time also constitutes a scenario where overestimation of the number of states is particularly prevalent.
Nevertheless, the ideas brought forward in this paper are not restricted to ecological applications, and are
potentially relevant in all fields of application of HMMs.
In HMM applications to animal movement data, the observed process is usually a bivariate time series
comprising step lengths and turning angles between subsequent locations, typically though not necessarily
collected using GPS technology. The states of the Markov chain underlying a fitted HMM can then be
interpreted as providing rough classifications of the behavioral modes of the animal observed (e.g. exploratory
6
vs. encamped; Morales et al., 2004; Patterson et al., 2016).
3 SIMULATION STUDIES
In this section, we present simulation studies to investigate the performance primarily of the AIC and the
BIC, but also of the ICL criterion, when it comes to selecting a suitable number of hidden states. Given the
asymptotic equivalence of AIC and cross-validated likelihood (Stone, 1977), we did not implement the latter
because of the associated substantial computational effort. ICL in contrast offers a conceptually different,
classification-based approach to assessing a model’s suitability, with a stronger focus on separation of classes,
such that it seems worthwhile to be investigated despite the fact that it is not regularly applied in practice.
Weshowcasescenarioswherethereisadditionalstructureinthedatathatisnotaccommodatedwithinbasic
HMM formulations. Each type of additional structure considered may be found in real data, and especially
within animal behavior data, where the assumptions made with the basic HMM formulation typically are
overly simplistic. What will be shown is that in such cases, the misspecification of the model formulation will
often be compensated by additional states which, to some extent, “mop up” the neglected structure. As a
consequence, standard model selection criteria (AIC and BIC) point to models with larger number of states
than necessary and adequate. The ICL criterion is much less prone to this behavior, but, as will be detailed
below, has other shortcomings.
We examine seven different scenarios: (1) outliers, (2) inadequate state-dependent distributions, (3) tem-
poralvariation,(4)individualheterogeneity,(5)asemi-Markovstateprocess,(6)asecond-orderMarkovchain
as state process, and (7) a violation of the conditional independence assumption. We start by describing the
scenarios in detail (Section 3.1), and subsequently provide an overview of the results and conclusions drawn
from these (Section 3.2).
3.1 SIMULATION SCENARIOS CONSIDERED
Asabaselinemodel,weconsideratwo-statestationarygamma–HMM,i.e.anHMMwithaunivariate,gamma-
distributedoutputvariable,withtwodistinctsetsofparameters,oneforeachstateoftheunderlyingstationary
Markov chain. This type of model is common for analyzing animal movement data where the observations
are the step lengths, with the two HMM states corresponding roughly to “foraging/resting” and “traveling”
behavior states, respectively. The two state-dependent gamma distributions in the baseline model are shown
in Figure 2. If not explicitly stated otherwise in the description of the individual scenarios below, then the
gammadistributionwithinthefirststatehasmean=0.5andshape=0.7(resultinginastrictlymonotonically
decreasingdensityfunction),andthegammadistributionwithinthesecondstatehasmean=4andshape=2.5
(resulting in a density function with mode distinct from zero).
In addition, if not explicitly stated otherwise, then the probability of leaving a state in any given time
7
6
0.
5
0.
4
0.
y state 1
ensit 0.3 state 2
d
2
0.
1
0.
0
0.
0 2 4 6 8 10
observation (e.g. step length)
Figure 2: Gamma distributions of the baseline HMM used in the simulation experiments.
interval is specified to be 0.1, leading to the transition probability matrix (t.p.m.)
 
0.9 0.1
Γ= .
0.1 0.9
In each of the seven scenarios considered below, some component of the baseline model formulation will
be slightly modified when simulating data. However, crucially, all the scenarios still involve only two genuine
states. We will then demonstrate that when not taking the modification into account, i.e. when fitting the
slightly misspecified basic models to the data, then additional states will be included as per recommendation
of model selection criteria, in order to compensate for the inflexibility of the model to otherwise capture the
additional structure.
Scenario 1 (outliers)
This first, very simple scenario represents a situation in which some of the data are outliers. In movement
ecology,theseoutlierscouldforexamplebeduetounusuallylargemeasurementerrors(e.g.asaresultofpoor
satellite coverage when using GPS tags). The corresponding simulated data are generated using the baseline
setup described above, but subsequently adding uniformly distributed random errors from the interval [10,20]
to only 0.5% of the data points ( 25 data points per sample generated). Intuitively it is clear that these few
≈
outlying values may cause the two-state baseline model to have a rather poor goodness of fit, since the two
gammastate-dependentdistributionsmaynotbeabletocovertheextremevalueswithoutlosingaccuracyfor
the non-outlying observations. This can potentially be compensated for by including additional states merely
to capture the outlying values.
Scenario 2 (inadequate emission distribution)
While parametric distributions will often provide good approximations of the actual empirical distribution
within a state, in practice it will practically never be the case that the true within-state distributions are
8
actual members of some parametric family. In other words, the parametric family being used (e.g. gamma
or normal distributions) will in general only provide an approximation of the actual data-generating process
within a state. In this scenario, we illustrate possible consequences of this for order selection.
To do so, the observations within the second state were not generated by a gamma distribution with mean
4andshape2.5anymore,butinsteadbyasimilarlyshapedbutheavy-taileddistributionwhichweconstructed
nonparametricallyusingB-splinebasisfunctions(Langrocket al.,2015). Thespecificshapeofthedistribution
within state 2 is shown in Figure 3.
0
3
0.
0
2
y 0.
sit
n
e
d 0
1
0.
0
0
0.
0 2 4 6 8 10 12
observation (e.g. step length)
Figure 3: Heavy-tailed emission distribution within state 2, as implemented in Scenario 2.
Clearly, the particular shape of this distribution cannot be fully captured by a single gamma distribution.
However,thedistributionappearstobesuchthatatwo-componentmixtureofgammadistributionswithinthe
second state may be sufficiently flexible to provide a good approximation to the nonparametric distribution.
(cid:0) (cid:1)
Notably though, a corresponding two-state model, with t.p.m. Γ = γ , i,j = 1,2, and mixture weights α
ij
and 1 α for the two gamma distributions within state 2, is equivalent to a three-state gamma–HMM with
−
suitably structured t.p.m.,
 
γ α(1 γ ) (1 α)(1 γ )
11 11 11
 − − − 
Γ=(1 γ22) αγ22 (1 α)γ22 .
 − − 
(1 γ ) αγ (1 α)γ
22 22 22
− −
Therefore, a two-state HMM with a gamma mixture in one of the states can just as well be represented by a
three-statesimplegamma–HMM.Itisthuspossiblethatmodelselectioncriteriafavormodelswithmorethan
two states, not because there are more than two genuine (behavioral) states, but because with the additional
states it is possible to represent more flexible emission distributions.
Scenario 3 (temporal variation)
In many ecological (and other) time series, there are clear temporal patterns in the data. We use diel (24 h
period) variation as an example but the issue applies to any temporal resolution with variation in the data
(e.g.seasonal, annual). Dielpatternscouldbepresentinthetransitionprobabilitiesand/orintheparameters
9
of the emission distributions, with the corresponding parameters then being cyclic functions of time.
The scenario we consider here is a hypothetical setting with a nocturnal animal that is more likely to be
active at nighttime than during the day. The state-switching probabilities will then depend on the time of
day. We constructed the transition probabilities using trigonometric functions, with a possible state switch
occurring every 15 minutes. The resulting transition probabilities, as a function of time of day, are shown in
Figure 4.
0
1. g12
g
21
s 8
bilitie 0.
a 6
ob 0.
pr
n 4
sitio 0.
tran 0.2
0
0.
0 4 8 12 16 20 24
time of day
Figure 4: transition probabilities as a function of time
In the movement ecology literature, analyses including such temporal components were presented for ex-
ample in Towner et al. (2016) and in Li and Bolker (in press). Although less clear as for example in Scenarios
2 and 3, it is possible that such a temporal structure in the data, if neglected in the model formulation, may
to some extent be captured by additional hidden states.
Scenario 4 (individual heterogeneity)
When observing more than only one individual, it is natural to assume the individuals to differ, inter alia,
in their personality and fitness, causing heterogeneity in the behavior between individuals. In the context of
animal movement, one could imagine younger individuals to move faster when in an active state, or to occupy
certainbehavioralstatesmoreoften,thanolderindividuals. Intheformercase,thenaturalmodelingapproach
forsuchasettingwouldbetoconsiderindividual-specificsetsofparametersofthegammadistributionwithin
the active state. For parsimony in terms of the number of model parameters, random effects would typically
be used (Schliehe-Diecks et al., 2012).
Not accounting for such individual heterogeneity within the model formulation, i.e. assuming identical
within-stategammadistributionsacrossallindividuals,couldagainleadtoinformationcriteriafavoringmodels
with more than two states: for example, one of the resulting >2 nominal HMM states may be capturing the
active movement of fitter individuals, while another may be associated with the active movement of less fit
individuals, with at least a third state capturing the inactive movement behavior. In our simulations, we
simulated 10 animal tracks of length 500 each, using a log-normal distribution with the parameters µ = ln4
10