Flexible Modelling of Diel and Other Periodic Variation in Hidden Markov Models

Animal behaviour is often characterised by periodic patterns such as seasonality or diel variation. Such periodic variation can be comprehensively studied from the increasingly detailed ecological time series that are nowadays collected, e.g. using GPS tracking. Within the class of hidden Markov models (HMMs), which is a popular tool for modelling time series driven by underlying behavioural modes, periodic variation is commonly modelled by including trigonometric functions in the linear predictors for the state transition probabilities. This parametric modelling can be too inflexible to capture complex periodic patterns, e.g. featuring multiple activity peaks per day. Here, we explore an alternative approach using penalised splines to model periodic variation in the state-switching dynamics of HMMs. The challenge of estimating the corresponding complex models is substantially reduced by the expectation–maximisation algorithm, which allows us to make use of the existing machinery (and software) for nonparametric regression. The practicality and potential usefulness of our approach is demonstrated in two real-data applications, modelling the movements of African elephants and of common fruit flies.


Introduction
Ecological time series data are often characterised by periodicities such as diel variation, i.e. recurrent patterns over a 24-h period. Ignoring periodic variation can invalidate statistical inference, e.g. standard errors might be underestimated due to residual autocorrelation [1]. Perhaps more importantly, adequately modelling such periodic variation is crucial to comprehensively understand behavioural dynamics, for example to identify times of day at which individuals tend to be most active, allowing inference on a species' temporal niche [2,3].
Fortunately, technological advances in, e.g. GPS tracking, accelerometry, and computer vision allow ecologists to study diel variation in much more detail than was previously possible [4,5]. One popular tool for modelling ecological time series data and the periodicities therein is given by the class of hidden Markov models (HMMs), which links the observed ecological data (e.g. step lengths and turning angles in animal movement) to underlying non-observable states (e.g. resting, foraging, travelling) [6].
In the existing literature, two different approaches have been used to infer periodic variation using HMMs. First, relatively basic HMMs can be used to infer an animal's behavioural sequence (state decoding), based on which diel variation can be investigated using simple visualisations [5,7,8] or an additional regression analysis [9,10]. From the statistical perspective, such a two-stage approach will often not be ideal: the uncertainty in the state allocation is not propagated, statistical inference on the periodic effects is not straightforward, and the dimension of the state space may be overestimated due to the misspecification of the basic model (see [11] for the latter). Second, periodic variation is nowadays often directly incorporated in HMMs using trigonometric modelling, for instance by relating the state transition probabilities to the hour of the day using sine and cosine functions with 24-h periods [12][13][14][15][16]. While this will often be sufficient, such a parametric modelling of the periodic effect may lack flexibility to capture complex periodic variation, for example with multiple activity peaks over the day. In principle, this limitation can be overcome by including multiple sine and cosine basis functions, with different wavelengths [17][18][19][20]. However, this can lead to numerical instability, and it can be tedious to select an adequate order.
In this contribution, we explore a more flexible, nonparametric estimation of periodicities in the state-switching dynamics of an HMM using cyclic splines. Thereby we avoid making any a priori assumptions on the functional shape of the periodic effect, thus allowing to infer arbitrarily complex behavioural diel patterns. For inference, we devise an expectation-maximisation (EM) algorithm, thereby isolating the estimation of the nonparametric periodic effect. This allows us to exploit the powerful machinery available for nonparametric (regression) modelling, specifically P-splines or other smoothing methods implemented in existing software packages such as mgcv [21]. The feasibility of the proposed approach is illustrated in two case studies, where we investigate diel variation of African elephants (Loxodonta africana) and of common fruit flies (Drosophila melanogaster).

Notation and Basics
HMMs are used to model time series data x 1 , . . . , x T (e.g. step lengths of an animal) driven by underlying states s 1 , . . . , s T (e.g. the behavioural modes). In a basic HMM, the latent state process is assumed to be a Markov chain with N states, characterised by the initial state probabilities δ (1) j = Pr(S 1 = j), j = 1, . . . , N , and the transition probability matrix (t.p.m.) Covariates-including time of day-can be included in either the state-dependent distributions f 1 , . . . , f N or the state transition probabilities γ (t) i j . We focus on the latter, as in ecological applications the main interest typically lies in the state process and its drivers, including temporal effects.

Trigonometric Modelling of Time-of-Day Variation
Including covariates in the state process amounts to regression modelling within the HMM, where covariates affect the state transition probabilities and hence, the behavioural decisions made by an animal. Specifically, if at time t − 1 the animal is in state i, the categorical distribution of states at time t is given by the vector (γ (t) i1 , . . . , γ (t) i N ). The covariate-dependence of this categorical distribution is typically modelled using a multinomial logistic regression, which is achieved by applying the inverse multinomial logit link to each row i of the t.p.m., defining τ (t) ii = 0 (reference category). The linear predictors τ (t) i j can include, inter alia, simple linear effects, polynomial effects, interaction terms, and random effects. Without the latter (for ease of notation), the general form of the linear predictor for γ (t) i j is where z tk are covariates. When the aim is to model periodic patterns in the stateswitching dynamics, the linear predictor can be extended by including trigonometric basis functions with the desired periodicity. For example, for modelling diel variation in a time series with hourly data, a possible simple form of the linear predictor is with the additional coefficients ω (i j) and ψ (i j) to be estimated alongside β (i j) . General periodicities are modelled analogously, replacing the 24 in the denominator by the period length (i.e. the number of sequential observations before completing one period). With only two harmonics, the flexibility of the periodic component of the linear predictor is somewhat limited: for example, this formulation implies that the periodic component has only one maximum turning point (such that patterns with multiple activity peaks throughout a day may not be adequately captured). The flexibility can be increased by including additional trigonometric functions with different wavelengths: see for example [19,20]. By increasing K , arbitrary (smooth) modelling of the periodic effect can be achieved. However, when complex periodic patterns are to be modelled, it can be tedious to select an adequate order K , with the risk of overfitting looming. It may then be more straightforward to avoid making any assumptions on the functional shape of the periodic effect, instead using nonparametric smoothing methods. Such an approach to estimating periodic effects will be presented and explored in the following.

Cyclic Splines for Modelling Time-of-Day Variation
We now consider nonparametric modelling of the periodic effect, replacing the sum of trigonometric basis functions in (2) by a spline. Specifically, we construct this spline as a linear combination of Q basis functions, with the scaling coefficients a Q to be estimated. We use cubic B-spline basis functions B 1 , ..., B Q , which are easy to compute and yield visually smooth functions. To enforce the desired periodicity, these are wrapped at the boundaries of the support [21]; see Fig. 1 for an illustration with Q = 8 and period length 24 h-again general periodicities are modelled analogously. In practice, a large Q (e.g. 20) is typically used to guarantee sufficient flexibility. Overfitting is avoided by including a penalty on the sums of squared differences between the coefficients a (i j) q associated with adjacent B-splines-an approach commonly referred to as P-spline modelling, cf. [22]. While strictly speaking, this model formulation is still parametric, it is commonly labelled as a nonparametric approach because with a large Q the modelling flexibility is effectively unlimited, and there is no meaningful interpretation of the coefficients a (i j) q .

EM-Based Estimation of HMMs with Cyclic Splines
The model formulation presented in the previous section effectively involves nonparametric regression modelling within HMMs. For example, in case of N = 2 states, the model features a nonparametric logistic regression for each of the state-switching probabilities γ (t) 12 and γ (t) 21 (see Eq. (1)). For such nonparametric regression modelling, the inferential machinery (including software packages) is well-established. Therefore, we apply the expectation-maximisation (EM) algorithm to isolate the estimation of the logistic regression component from the estimation of the other parameters of the HMM, in particular those associated with the state-dependent process. This allows us to exploit the tools available for nonparametric logistic regression modelling.
To set up the EM algorithm, we consider the complete-data likelihood (CDLL) of the HMM, i.e. the joint log-likelihood of the observations and the states, with θ the set of parameters necessary to define δ (1) j , (t) , and the state-dependent distributions f j (x). Each iteration of the EM algorithm involves an E-step, replacing all functions of the unobserved states in the CDLL by their conditional expectations (given the data and the current parameter values), and an M-step, optimising the resulting CDLL with respect to θ . In the case of HMMs, the appeal of the EM algorithm lies in the fact that the M-step neatly splits into several separate optimisation problems-namely one for each the initial distribution, the t.p.m., and the state-dependent distributionswhich we exploit below to conveniently estimate the cyclic spline component.
To apply the E-step, we define the indicator variables and rewrite the CDLL as In the E-step, the indicator variables are then replaced by their conditional expectationŝ with θ the current guess of the parameter vector. These conditional expectations are calculated using the standard forward and backward recursions (see [23]). The M-step then involves optimising the CDLL (5), with u i (t) and v i j (t) replaced byû i (t) andv i j (t), respectively, with respect to the parameter vector θ. The updated estimates of the initial state distribution as well as the state-dependent distribution are obtained as comprehensively described in [23]. For the particular model formulation considered in this contribution, the interest (and challenge) lies in the update of the parameters that affect the state transition probabilities, i.e. the second term in (5), The first summation in (6) corresponds to the N rows of the t.p.m., each of which implies a categorical regression model for the transition to the next state. We can estimate the associated parameters of each of these regressions separately. For example, Each of these two terms is the log-likelihood of a logistic regression model. For example, the first term can be rewritten as In logistic regression terminology, the sum T t=2 v 12 (t) gives the number of "successes" (here, the number of switches from state 1 to state 2), and T t=2 v 11 (t) is the number of "failures" (here the number of instances when the process remains in state 1). Within EM, the indicator variables v i j (t) are replaced by their conditional expectationsv i j (t), such that (7) becomes a weighted log-likelihood.
The time-varying transition probability γ (t) 12 in (7) is modelled using cyclic Psplines, see (1) and (3). As described in Sect. 2.3, a wiggliness penalty is added to the weighted log-likelihood, for example with 2 denoting the second-order difference operator and λ i the (state-dependent) smoothing penalty (see [22]). The estimation of this weighted nonparametric logistic regression can conveniently be conducted using well-established machinery, including existing software. In the case studies below, we implemented this part of the M-step in the EM algorithm using the mgcv package in R [21]. The updated parameter estimates are then used in the E-step of the next iteration. The E and M steps are repeated until a convergence criterion defined by the user is reached [24], e.g. that the difference between the likelihood values obtained in two consecutive iterations is below some threshold. This iterative scheme identifies a (local) maximum of the likelihood function. To increase the chances of finding the global maximum, several initial starting values for θ should be tested.

African Elephant
We consider hourly GPS data collected for an African elephant in Etosha National Park, Namibia, from October 2008 to August 2010. The data are available from the Movement Bank Repository [25,26], cf. [27]. From the positional data, we calculate the Euclidean step lengths as well as the turning angles between consecutive compass directions. From these two metrics, we aim to investigate diel patterns in the elephant's behaviour. The empirical step length distributions throughout the day indicate a relatively complex diel variation with two activity modes, which may be difficult to adequately model using parametric periodic effects (Fig. 2).
We model the data using 2-state HMMs with gamma and von Mises distributions for the step lengths and turning angles, respectively, assuming conditional independence of the two variables, given the states [28,29]. For modelling diel variation in the state-switching dynamics, we consider the cyclic P-spline approach (using the default options implemented in mgcv), the trigonometric approach (2) with K = 1, 2 and (3), and, as an additional benchmark, a model with homogeneous Markov chain (i.e. no diel variation). All fitted models feature an "encamped" state with relatively short step lengths and frequent reversals in direction (state 1) and an "exploratory" state with longer steps and higher persistence in direction (state 2)-see Fig. 6 in Appendix. Figure 3 displays the time-varying probabilities of switching states (left panel: from state 1 to state 2, right panel: vice versa) as estimated under the nonparametric as well as the parametric approach. All models detect a reduction in exploratory activity during the night. However, the flexible P-spline approach additionally captures a bimodal diel variation, with more frequent switching to the exploratory mode in the early morning hours but also in the early afternoon. In contrast, the commonly used trigonometric effect modelling with K = 1 (i.e. one sine and one cosine basis function) is not sufficiently flexible to identify this bimodality. When increasing the order to K = 2, the bimodality can be identified, however, only with K = 3 the parametric approach produces results similar to those obtained using splines (thus indicating that even with K = 2 the parametric model might be too inflexible). Furthermore, the proportion of time spent in the exploratory state, for each time of day calculated based on the Viterbi-decoded states, varies notably across the five models fitted (see Fig. 4). This underlines the importance of adequately modelling diel variation, as inflexible models can invalidate inference on the state process.  To formally compare spline-based models with parametric alternatives using trigonometric base functions, we consult the Akaike information criterion (AIC) and the Bayesian information criterion (BIC; see Table 1). The AIC and the BIC favour the trigonometric models with K = 5 and K = 3, respectively, with the spline-based approach arriving at a similarly complex model as measured by the effective degrees of freedom (edf). Notably, this demonstrates that the penalised spline approach corresponds to data-driven model selection since the choice of the wiggliness penalty λ aims at achieving a favourable balance between underfitting and overfitting.

Common Fruit Flies
In the second case study, we consider the locomotor activity of laboratory wild type Drosophila melanogaster (iso31) [30]. We collected 2-to 3-days-old male flies and trained them individually to a standard 12-h-light and 12-h-dark condition (LD) for 4.5 days in locomotor tubes. Subsequently, we subjected them to 6 days of constant darkness (DD). The temperature was kept constant (25 • C). During these 10 days, locomotor activity was recorded using the Drosophila Activity Monitor (DAM) system (TriKinetics Inc), by counting the times each fly interrupts the infrared beam passing the middle of the locomotor tube. We consider two time series-one under light condition LD and the other under condition DD-for each of 15 individuals. Each observation is the count of beams crossed over a period of 30 min. The time series of half-hourly counts are modelled using a 2-state HMM with negative binomial state-dependent distributions. For the state transition probabilities, we consider the same time-varying predictors as in the elephant example, additionally allowing for different periodic effects under the two light conditions. The fitted models' states are associated with low and high activity, with state-dependent mean counts of 2.7 and 54.9, respectively, obtained for the spline-based model (see Fig. 7 in Appendix.). Figure 5 shows the time-varying probability of occupying the high-activity state as obtained for the different models considered. Note that these are effectively summary statistics implied by the time-varying t.p.m., which are shown here to facilitate the comparison between the two light conditions. The results emphasise the importance of allowing for sufficient modelling flexibility in the periodic effects, as the commonly used approach with only one sine and one cosine basis function fails to capture several key characteristics: (1) the bimodal activity pattern over the course of the day; (2) the near-certain occupancy of the high-activity state in the evening hours (and of the lowactivity state during the night) in the DD condition; (3) the fact that the first activity peak is less pronounced, and the second more pronounced, in the DD condition. The other three models-i.e. those with trigonometric effect modelling and either K = 2 or K = 3 as well as the spline-based model-all yield similar results. The slight midday  Models favoured by the criteria are highlighted in bold font peak only revealed by the spline-based approach-as well as by trigonometric models with K ≥ 4 (not displayed in the figure)-was also found in another study on activity patterns of fruit flies, albeit under varying temperature conditions [31]. Furthermore, the model comparison shows that the spline-based approach leads to a similar fit as the trigonometric approach, with the BIC favouring the latter with K = 4 (see Table 2).
patterns, e.g. multiple activity peaks throughout the day, we explored a flexible nonparametric approach using cyclic P-splines. As illustrated in two case studies, such flexibility in the modelling of periodic variation can uncover relevant patterns that may otherwise go unnoticed, emphasising the potential usefulness of the approach, in particular in settings where periodic variation is of primary interest (cf. [32][33][34]). We implemented the HMM with cyclic P-splines building on the mgcv functionality within the EM algorithm. However, this is not the only option for conducting inference for such a model. In particular, optimisation of the HMM's marginal likelihood, obtained by integrating out the spline coefficients using the Laplace approximation [35], has recently been explored and is already implemented in the immensely flexible R package hmmTMB [36]. Moreover, both hmmTMB and the presented EM algorithm are not limited to P-spline modelling but could also incorporate other functional relationships within HMMs, such as random effects or more complex smoothing functions provided by mgcv.
The manifold possibilities of flexibly modelling covariate effects further complicate the already difficult task of model selection in HMMs, which is aggravated by the interplay of the number of states and the modelling of the state process [37]. Although information criteria offer some guidance in identifying an adequately complex model, additional considerations regarding interpretability and computational costs need to be taken into account. For example, even if a spline-based model fits the data slightly better than a simpler parametric one, this advantage may be outweighed by the additional computational burden. Therefore, HMMs demand a holistic approach to model selection, that is to pragmatically balance goodness of fit and study aims.
In practice, nonparametric modelling of periodic variation allows to investigate temporal patterns without making any restrictive assumptions a priori. When used as an exploratory tool, the approach may of course also show that simple trigonometric modelling is sufficient. In both case studies presented in this contribution, information criteria did indeed favour trigonometric modelling of the periodic variation, though only when using considerably more basis functions than are commonly applied in ecological modelling. Selecting an adequate number of basis functions can be tedious in practice, such that the spline-based approach, with its automated optimisation of the bias-variance trade-off via penalised likelihood, may sometimes be more convenient to implement. In any case, incorporating the required flexibility to capture periodic variation-whether using nonparametric or parametric modelling-is crucial for guaranteeing valid inference on behavioural processes.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Appendix A
See Figs. 6 and 7.