\documentclass[12pt]{article}
\usepackage[latin9]{inputenc}
\usepackage{geometry}
\geometry{verbose}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{setspace}
\usepackage{natbib}
\setcitestyle{aysep={}}
\usepackage{adjustbox}
\usepackage[labelfont=bf,labelsep=period,font={small,sl}]{caption}
\usepackage{amsthm}
\usepackage{epsfig}
\usepackage{psfrag}
\usepackage{booktabs}
\usepackage{url}

% \usepackage{lineno}

\renewcommand{\baselinestretch}{1.3}

\begin{document}

<<include=FALSE>>=
library(knitr)
opts_chunk$set(
concordance=TRUE
)
@

%\VignetteIndexEntry{Guide to using momentuHMM}
%\VignetteEngine{knitr::knitr}

% set margin to 4cm for title page
\newgeometry{margin=4cm}

\begin{center}
  \texttt{\LARGE momentuHMM}{\LARGE : R package for analysis of telemetry data using generalized multivariate hidden Markov models of animal movement}\vspace{0.5in}
  \par
\end{center}

\begin{center}
  {\large Brett T. McClintock$^{1}$ and Th\'eo Michelot$^{2}$} 
  \par
\end{center}

\begin{center}
  \hrulefill{} 
  \par
\end{center}

\begin{center}
  \global\long\def\baselinestretch{1.25}
  {\large $^{1}$Marine Mammal Laboratory}\\
  {\large Alaska Fisheries Science Center}\\
  {\large {} NOAA National Marine Fisheries Service}\\
  {\large {} Seattle, U.S.A.}\\
  {\large {} {\em Email:} brett.mcclintock@noaa.gov} 
  \par
\end{center}

{\large \par}

\begin{center}
  {\large $^{2}$School of Mathematics and Statistics}\\
  {\large University of Sheffield}\\
  {\large {} Sheffield, U.K.}
  \par
\end{center}

\begin{center}
  {\large \hrulefill{}} 
  \par
\end{center}

\begin{center}
  \textsc{Running Head}: R package \verb|momentuHMM| \bigskip{}  
  \par
\end{center}

\begin{center}
  \today
  \par
\end{center}

\clearpage{}

% \setlength{\textheight}{575pt} \global\long\def\baselinestretch{2}
%% ABSTRACT %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%  make sure that the document has 25 lines per page (it is 12 pt)
% \setlength{\textheight}{575pt} \setlength{\baselineskip}{24pt} %think this may do doublespacing (required for submission I think)

\newpage{}

%\linenumbers

% set margin to 3cm for main text
\newgeometry{margin=3cm}

\noindent \textbf{Summary}\\
\textbf{1.} Discrete-time hidden Markov models (HMMs) have become an immensely popular tool for inferring latent animal behaviors from telemetry data, largely because they are relatively fast and easy to implement when data streams are observed without error and at regular time intervals. While movement HMMs typically rely solely on location data, auxiliary biotelemetry and environmental data are powerful and readily-available resources for incorporating much more behavioral realism and inferring ecological relationships that would otherwise be difficult or impossible to infer from location data alone.  However, there is a paucity of generalized user-friendly software available for implementing (multivariate) HMMs of animal movement. Furthermore, location measurement error, temporal irregularity, and other forms of missing data are often pervasive in telemetry studies (particularly in marine systems).\\ %and the incorporation of uncertainty attributable to location measurement error or missing data typically requires fitting HMMs using custom and computationally-demanding model fitting techniques, such as Markov chain Monte Carlo.\\ 
\textbf{2.} Here we provide a guide to using an open-source R package, \verb|momentuHMM| version \Sexpr{installed.packages()["momentuHMM","Version"]}%(Maximum likelihood analysis Of animal MovemENT behavior Using multivariate Hidden Markov Models)
, that addresses many of the deficiencies in existing software.  Features for multivariate HMMs in \verb|momentuHMM| (pronounced ``momentum'') include: 1) tools for data pre-processing and visualization; 2) user-specified probability distributions for an unlimited number of data streams and latent behavior states, such as those based on location (e.g., step length, turning angle) and auxiliary biotelemetry data (e.g., from pressure, conductivity, heart rate, or motion sensors); 3) biased and correlated random walk movement models, including ``activity centers'' associated with attractive or repulsive forces; 4) user-specified design matrices and constraints for covariate modelling of initial distribution, state transition probability, and probability distribution parameters using linear model formulas familiar to most R users; 5) multiple imputation methods that account for observation error attributable to measurement error and temporally-irregular or missing data; 6) seamless integration of spatio-temporal covariate raster data; 7) cosinor and spline regression formulas for cyclical (e.g., daily, seasonal) and other complicated patterns; 8) model checking and selection; and 9) data simulation capabilities for study design, power analyses and assessing model performance, including simulation of location data subject to temporal irregularity and/or measurement error.\\ 
\textbf{3.} After providing a brief introduction to (multivariate) HMMs for telemetry data, we demonstrate some of the capabilities of \verb|momentuHMM| using real-world examples. This brief tutorial includes workflows for data formatting, model specification, model fitting, and diagnostics.\\
\textbf{4.} While many of the features of \verb|momentuHMM| were motivated by animal movement data, the package can be used for analyzing any type of data that is amenable to (multivariate) HMMs. Practitioners interested in additional features for \verb|momentuHMM| are encouraged to contact the authors.\\


\noindent \textbf{Key-words} animal biotelemetry, biologging, \verb|crawl|, \verb|moveHMM|, state-space model, state-switching

%\newpage{}

\global\long\def\baselinestretch{1.0}
 \global\long\def\baselinestretch{1.0}


\section{Introduction}

Discrete-time hidden Markov models (HMMs) have become immensely popular for the analysis of animal telemetry data \citep[e.g.][]{MoralesEtAl2004,JonsenEtAl2005,LangrockEtAl2012,McClintockEtAl2012}. In short, an HMM is a time series model composed of a (possibly multivariate) observation process $({\mathbf Z}_1,\ldots,{\mathbf Z}_T)$, in which each data stream is generated by $N$ state-dependent probability distributions, and where the unobservable (hidden) state sequence $(S_t\in\{1,\ldots,N\},t=1,\ldots,T)$ is assumed to be a Markov chain.  The state sequence of the Markov chain is governed by (typically first-order) state transition probabilities, $\gamma_{ij}^{(t)}=\text{Pr}(S_{t+1}=j \mid S_t=i)$ for $i,j=1,\ldots,N$, and an initial distribution ${\boldsymbol \delta}^{(0)}$.  The likelihood of an HMM can be succinctly expressed using the forward algorithm:
\begin{equation}
  \mathcal{L}={\boldsymbol \delta}^{(0)} {\mathbf \Gamma}^{(1)} {\mathbf P}({\mathbf z}_1) {\mathbf \Gamma}^{(2)} {\mathbf P}({\mathbf z}_2) {\mathbf \Gamma}^{(3)} \cdots {\mathbf \Gamma}^{(T-1)} {\mathbf P}({\mathbf z}_{T-1}) {\mathbf \Gamma}^{(T)} {\mathbf P}({\mathbf z}_{T}) {\bf 1}^N,
  \label{eq:HMMlike}
\end{equation}
where ${\mathbf \Gamma}^{(t)}=\left(\gamma_{ij}^{(t)} \right)$ is the $N \times N$ transition probability matrix, ${\mathbf P}({\mathbf z}_t)=\text{diag}(p_1({\mathbf z}_t), \ldots, p_N({\mathbf z}_t))$, $p_s({\mathbf z}_t)$ is the conditional probability density of ${\mathbf Z}_t$ given $S_t=s$, and ${\bf 1}^N$ is a $N$-vector of ones \citep[for a thorough introduction to HMMs see][]{ZucchiniEtAl2016}.  

One of the most common discrete-time animal movement HMMs for telemetry location data is composed of two data streams, step length and turning angle (or bearing), which are calculated for each of the $T$ time steps from the temporally-regular observations of an animal's position, $(x_t,y_t)$, for $t=1,\ldots,T+1$ \citep[e.g.][]{MoralesEtAl2004,LangrockEtAl2012,McClintockEtAl2012}. Step length $(l_t)$ is typically calculated as the Euclidean distance between the locations $(x_t,y_t)$ and $(x_{t+1},y_{t+1})$, while turning angle $(\phi_t)$ is calculated as the change in bearing $\left(b_t=\text{atan2}(y_{t+1}-y_t,x_{t+1}-x_t)\right)$ between the intervals $[t-1,t]$ and $[t,t+1]$ (e.g. $\phi_t=0$ if $b_{t-1}=b_t$). For this HMM composed of 2 data streams, ${\mathbf z}_t=(l_t,\phi_t)$, and, conditional on the latent state $S_t$, independent probability distributions are typically assumed for each stream; that is, $p_s({\mathbf z}_t)=p_s(l_t)p_s(\phi_t)$. Some common probability distributions for the step length data stream are the gamma or Weibull distributions, while the wrapped Cauchy or von Mises distributions are often employed for turning angle or bearing. For a fitted HMM, the Viterbi algorithm is used to compute the most likely sequence of underlying states \citep{ZucchiniEtAl2016}. In movement HMMs, the states are often considered as proxies for animal behaviour.

While HMMs for animal movement based solely on location data are somewhat limited in the number and type of biologically-meaningful movement behavior states they are able to accurately identify, advances in biologging technology are now allowing the collection of valuable auxiliary biotelemetry data (e.g., dive activity, accelerometer, heart rate, stomach temperature), which, when combined with location data, allow for multivariate HMMs that can incorporate much more behavioral realism and facilitate inferences about complex ecological relationships that would otherwise be difficult or impossible to infer from location data alone \citep[e.g.][]{McClintockEtAl2013c,DeRuiterEtAl2017,McClintockEtAl2017}.  Multivariate HMMs that utilize both location and auxiliary biotelemetry data can facilitate the identification of additional states that go beyond the $N=2$ state approaches that are most frequently used by practitioners. For example, the most widely used 2-state HMMs for animal movement include ``encamped'' (or ``foraging'') and ``exploratory'' (or ``transit'') states characterized by area-restricted-search-type movements (shorter step lengths with little to no directional persistence) and migratory-type movements (longer step lengths with high directional persistence), respectively \citep{MoralesEtAl2004,JonsenEtAl2005}.  However, very different behaviors can exhibit similar horizontal trajectories.  For example, for herbivores such as North American elk \citep{MoralesEtAl2004} or central-place foragers such as harbour seals \citep{McClintockEtAl2013c}, the horizontal trajectories of ``resting'' and ``foraging'' movements can be very difficult to distinguish. Standard 2-state HMMs based solely on horizontal trajectory will tend to lump these behaviors together, and this could have unintended consequences if, for example, one intends to use the estimated state sequences to identify foraging habitat. In order to tweeze apart distinct behaviors with similar horizontal trajectories, additional states can be informed by auxiliary information (such as mandible accelerometer or dive data), incorporated as additional data stream(s) in a multivariate HMM.

When data streams are observed without error and at regular time intervals, a major advantage of HMMs is the relatively fast and efficient maximization of the likelihood using the forward algorithm (Eq. \ref{eq:HMMlike}).  However, location measurement error is rarely non-existent in animal-borne telemetry studies and depends on both the device and the system under study.  For example, GPS errors are typically less than 50m, but Argos errors can exceed 10km \citep[e.g.][]{CostaEtAl2010}.  An extreme case of missing data can arise when location data are obtained with little or no temporal regularity, as in many marine mammal telemetry studies \citep[e.g.][]{JonsenEtAl2005}, such that few (if any) observations align with the regular time steps required by discrete-time HMMs. When explicitly accounting for uncertainty attributable to location measurement error, temporally-irregular observations, or other forms of missing data, one must typically fit (multivariate) HMMs using computationally-intensive (and often time-consuming) model fitting techniques such as Markov chain Monte Carlo \citep{JonsenEtAl2005,McClintockEtAl2012}.  However, complex analyses requiring novel statistical methods and custom model-fitting algorithms are not practical for many practitioners.

While statisticians have been applying HMMs to telemetry data for decades, R \citep{RCoreTeam2017} packages such as \verb|bsam| \citep{JonsenEtAl2005}, \verb|moveHMM| \citep{MichelotEtAl2016}, and \verb|swim| \citep{WhoriskeyEtAl2017} have recently helped make these models of animal movement behavior more accessible to the practitioners that are actually conducting telemetry studies.  These advances represent important steps toward making HMMs of animal movement more accessible, but the models that can currently be implemented using existing software remain limited in many key respects. For example, existing HMM software for animal movement is limited to two data streams based solely on location data (e.g. step length and turning angle), and while \verb|moveHMM| allows for a user-specified number of latent behavioral states (\verb|bsam| and \verb|swim| are limited to $N=2$ states), it is typically difficult to identify $>$2 biologically-meaningful behavior states from only 2 data streams \citep[e.g.][]{MoralesEtAl2004,BeyerEtAl2013,McClintockEtAl2014b}. Both \verb|moveHMM| and \verb|swim| are designed for temporally-regular (or linearly-interpolated) location data with negligible measurement error, but the realities of animal-borne telemetry often yield temporally-irregular location data subject to error (particularly in aquatic environments). Other notable deficiencies of existing software include limited abilities to incorporate spatio-temporal environmental or individual covariates on parameters, biased (or directed) movements in response to attractive or repulsive forces \citep[e.g.][]{McClintockEtAl2012,LangrockEtAl2014}, cyclical (e.g. daily, seasonal) and other more complicated behavioral patterns, or constraints on parameters. 

To address these deficiencies in existing software, we developed a user-friendly R package, \verb|momentuHMM| (Maximum likelihood analysis Of animal MovemENT behavior Using multivariate Hidden Markov Models), intended for practitioners wishing to implement more flexible and realistic (multivariate) HMM analyses of animal movement while accounting for common challenges associated with telemetry data \citep{McClintockMichelot2018}. Features for multivariate HMM analyses in \verb|momentuHMM| include: 1) tools for data pre-processing and visualization; 2) user-specified probability distributions for an unlimited number of data streams and latent behavior states; 3) biased and correlated random walk movement models, including ``activity centers'' associated with attractive or repulsive forces \citep[e.g.][]{McClintockEtAl2012}; 4) user-specified design matrices and constraints for covariate modelling of state transition probability and probability distribution parameters using linear model formulas familiar to most R users; 5) multiple imputation methods that account for observation error attributable to measurement error and temporally-irregular or missing data \citep{HootenEtAl2017,McClintock2017}; 6) seamless integration of spatio-temporal environmental covariate data (e.g., wind direction, forest cover, sea ice concentration) using the \verb|raster| package \citep{Hijmans2016}; 7) cosinor \citep[e.g.][]{Cornelissen2014} and spline regression formulas for cyclical and other complicated behavioral patterns; 8) model checking and selection; and 9) data simulation capabilities for study design, power analyses and assessing model performance, including simulation of location data subject to temporal irregularity and/or measurement error. 

In the following tutorial, we demonstrate some of the capabilities of \verb|momentuHMM| using real-world examples, including an example of periodic cycles in African elephant movement, a 3-state (``resting'', ``foraging'', ``transit'') northern fur seal example incorporating auxiliary dive activity data \citep{McClintockEtAl2014b}, a loggerhead turtle example relating ``foraging'' and ``transit'' movements to ocean surface currents, a 5-state grey seal example incorporating biased movements toward haul-out and foraging locations \citep{McClintockEtAl2012}, a 4-state (``outbound'', ``searching'', ``foraging'', ``inbound'') southern elephant seal example with biased movements toward and away from a colony \citep{MichelotEtAl2017}, a 3-state (``resting'', ``foraging'', ``transit'') harbour seal example using population-level constraints on movement parameters \citep{McClintockEtAl2013c}, and a 6-state northern fulmar example incorporating biased movements relative to both static (i.e. colony) and dynamic (i.e. fishing vessels) activity centers \citep{PirottaEtAl2018}. Using simulated data, we also demonstrate how the group dynamic model of \cite{LangrockEtAl2014} can be implemented using \verb|momentuHMM|. This brief tutorial includes workflows for data formatting, model specification, model fitting, and diagnostics. While many of the features of \verb|momentuHMM| were motivated by animal movement data, the package can be used for analyzing any type of data that is amenable to (multivariate) HMMs.  Additional information, including help files, data, examples, and package usage is available by downloading the \verb|momentuHMM| package from CRAN (\url{cran.r-project.org}) or Github (\url{github.com/bmcclintock/momentuHMM}). This article describes \verb|momentuHMM| version \Sexpr{installed.packages()["momentuHMM","Version"]}.

\section{momentuHMM overview}
%Given the deliberately similar syntax between the two packages, 
Before delving into some of the finer details, we will first provide an overview of the main features and functions of \verb|momentuHMM| (pronounced ``momentum''). While space is limited in this tutorial, further details on implementation can be found in the package's documentation and vignette. The workhorse functions of \verb|momentuHMM| are listed in Table \ref{tab:functions}. Usage of several of these functions (e.g. \verb|fitHMM|, \verb|prepData|, \verb|simData|) is deliberately very similar to equivalent functions in \verb|moveHMM| \citep{MichelotEtAl2016} , but the \verb|momentuHMM| arguments for these functions have been generalized and expanded to accommodate a more flexible framework for data pre-processing, model specification, parameterization, and simulation. R users already familiar with \verb|moveHMM| will therefore likely find it easy to immediately begin using \verb|momentuHMM|. %While the \verb|momentuHMM| syntax is therefore more complicated, any \verb|moveHMM| model can be implemented in \verb|momentuHMM|.
\begin{table}
  \caption{\label{tab:functions} Workhorse functions for the R package momentuHMM.}
  \begin{tabular}{ll}
  \toprule
  Function & Description \tabularnewline
  \midrule
  %\verb|AIC.momentuHMM| & AIC for one or several \verb|momentuHMM| models  \tabularnewline
  %\verb|CIbeta| & Confidence intervals for working (beta) parameters  \tabularnewline
  %\verb|CIreal| & Confidence intervals for natural (real) parameters  \tabularnewline
  \verb|crawlMerge| & Merge \verb|crawlWrap| output with additional data streams or covariates  \tabularnewline 
  \verb|crawlWrap| & Fit \verb|crawl| models and predict temporally-regular locations  \tabularnewline  
  \verb|fitHMM| & Fit a (multivariate) HMM to the data  \tabularnewline  
  \verb|MIfitHMM| & Fit (multivariate) HMMs to multiple imputation data  \tabularnewline  
  \verb|MIpool| & Pool \verb|momentuHMM| model results across multiple imputations  \tabularnewline 
  \verb|plot.crwData| & Plot \verb|crawlWrap| output \tabularnewline 
  \verb|plot.miSum| & Plot summaries of multiple imputation \verb|momentuHMM| models  \tabularnewline 
  \verb|plot.momentuHMM| & Plot summaries of \verb|momentuHMM| models  \tabularnewline 
  \verb|plot.momentuHMMData| & Plot summaries of selected data streams and covariates  \tabularnewline 
  \verb|plotPR| & Plot time series, qq-plots and sample ACFs of pseudo-residuals \tabularnewline 
  \verb|plotSat| & Plot locations on satellite image \tabularnewline   
  \verb|plotSpatialCov| & Plot locations on raster image \tabularnewline   
  \verb|plotStates| & Plot the (Viterbi-)decoded states and state probabilities \tabularnewline 
  \verb|plotStationary| & Plot stationary state probabilities \tabularnewline
  \verb|prepData| & Pre-process data streams and covariates \tabularnewline 
  \verb|pseudoRes| & Calculate pseudo-residuals for \verb|momentuHMM| models \tabularnewline 
  \verb|simData| & Simulate data from a (multivariate) HMM \tabularnewline 
  \verb|stateProbs| & State probabilities for each time step \tabularnewline 
  \verb|viterbi| & Most likely state sequence (using the Viterbi algorithm)  \tabularnewline  
  \bottomrule
  \end{tabular}
\end{table}

One of the key features of \verb|momentuHMM| is the ability to include an unlimited number of HMM data streams (e.g.\ step length, turning angle, dive activity, heart rate) arising from a broad range of commonly used probability distributions (e.g.\ beta, gamma, normal, Poisson, von Mises, Weibull). Any of the parameters of the probability distributions used for the observed data can be modelled as a function of environmental and individual covariates using link functions (Table \ref{tab:pdfs}). For any given ``natural scale'' (or ``real scale'') probability distribution parameter $\theta$, all of the link functions $(g)$ in \verb|momentuHMM| are of the general form $g({\boldsymbol \theta}) =  {\mathbf X}_\theta{\boldsymbol \beta}_\theta$, where ${\mathbf X}_\theta$ is the $T \times K$ design matrix (composed of $K$ covariates) and ${\boldsymbol \beta}_\theta$ is the correponding $K$-vector of ``working scale'' (or ``beta scale'') parameters for $\theta$. For example, suppose step length is assumed to have a gamma distribution, $l_t\mid S_t=s \sim \text{gamma}(\mu_s,\sigma_s)$. In \verb|momentuHMM|, the natural scale parameters for the gamma distribution are the (state-dependent) step length mean $(\mu_s>0)$ and standard deviation $(\sigma_s>0)$.  Because both of these parameters must be positive, the log link function is a natural choice for modelling these parameters as a function of covariates, e.g., $\log({\boldsymbol \mu}) =  {\mathbf X}_\mu  {\boldsymbol \beta}_\mu$ and $\log({\boldsymbol \sigma}) =  {\mathbf X}_\sigma  {\boldsymbol \beta}_\sigma$.

The state transition probabilities $({\mathbf \Gamma}^{(t)})$ and initial distribution $({\boldsymbol \delta}^{(0)})$ can also be modelled as functions of covariates, using a multinomial logit link, as described e.g.\ by \cite{MichelotEtAl2016}. Permissable R classes for covariates include \verb|numeric|, \verb|integer|, or \verb|factor|. Factors can be particularly useful for specifying models with individual- or group-level (e.g. sex or age class) effects on state transition and probability distribution parameters. Spatio-temporal covariates can also be of classes \verb|rasterLayer|, \verb|rasterStack|, or \verb|rasterBrick| \citep{Hijmans2016}, in which case \verb|momentuHMM| automatically extracts the appropriate covariate values from the raster based on the time and location of each observation (see example in section \ref{sec:turtle}).

\begin{table}
  \caption{\label{tab:pdfs} Data stream $(z)$ probability distributions, natural parameters, and default link functions for covariate modelling. Probability distributions with positive support can be zero-inflated (with additional zero-mass parameters), while the beta distribution can be zero- and/or one-inflated (with additional one-mass parameters). If user-specified bounds are provided, then custom link functions are used instead of the defaults (see package documentation for further details). If both zero- and one-inflation are included, then a multinomial logit link is used because these probabilities must sum to less than one (in this case, any user-specified bounds for the zero- and one-inflation parameters are ignored). If circular-circular regression is specified for the mean of angular distributions (``vm'' and ``wrpcauchy''), then a link function based on \cite{RivestEtAl2016} is used. The von Mises consensus distribution (``vmConsensus'') is a von Mises circular-circular regression model where the concentration parameter depends on the level of agreement among short-term directional persistence and angular covariates. Users seeking additional probability distributions are encouraged to contact the authors.}
  \begin{tabular}{llll}
  \toprule
  Distribution & Support & Parameters & Link function\footnotemark \tabularnewline
  \midrule
  Bernoulli (``\verb|bern|'')                 & $z_t\in\{0,1\}$          & $\verb|prob|\in(0,1)$            &  $\text{logit}$ \tabularnewline  
  Beta (``\verb|beta|'')                      & $z_t\in(0,1)$            & $\verb|shape1|>0$                &  $\log$ \tabularnewline  
                                              &                          & $\verb|shape2|>0$                &  $\log$ \tabularnewline
                                              &                          & $\verb|zero-mass|\in(0,1)$       &  $\text{logit}$ \tabularnewline 
                                              &                          & $\verb|one-mass|\in(0,1)$        &  $\text{logit}$ \tabularnewline 
  Exponential (``\verb|exp|'')                & $z_t>0$                  & $\verb|rate|>0$                  &  $\log$ \tabularnewline  
                                              &                          & $\verb|zero-mass|\in(0,1)$       &  $\text{logit}$ \tabularnewline 
  Gamma (``\verb|gamma|'')                    & $z_t>0$                  & $\verb|mean|>0$                  &  $\log$ \tabularnewline  
                                              &                          & $\verb|sd|>0$                    &  $\log$ \tabularnewline  
                                              &                          & $\verb|zero-mass|\in(0,1)$       &  $\text{logit}$ \tabularnewline 
  Log normal (``\verb|lnorm|'')               & $z_t>0$                  & $\verb|location|\in{\rm I\!R}$   &  identity \tabularnewline  
                                              &                          & $\verb|scale|>0$                 &  $\log$ \tabularnewline  
                                              &                          & $\verb|zero-mass|\in(0,1)$       &  $\text{logit}$ \tabularnewline 
  Normal (``\verb|norm|'')                    & $z_t\in{\rm I\!R}$       & $\verb|mean|\in{\rm I\!R}$       &  identity \tabularnewline  
                                              &                          & $\verb|sd|>0$                    &  $\log$ \tabularnewline 
  Poisson (``\verb|pois|'')                   & $z_t\in\{0,1,\ldots\}$   & $\verb|lambda|>0$                &  $\log$ \tabularnewline  
  Von Mises (``\verb|vm|'')                   & $z_t\in(-\pi,\pi]$       & $\verb|mean|\in(-\pi,\pi]$       &  $\tan(\verb|mean|/2)$ \tabularnewline  
                                              &                          & $\verb|concentration|>0$         &  $\log$ \tabularnewline 
  Von Mises (``\verb|vmConsensus|'')          & $z_t\in(-\pi,\pi]$       & $\verb|mean|\in(-\pi,\pi]$       &  \cite{RivestEtAl2016} \tabularnewline  
                                              &                          & $\verb|kappa|>0$                 &  $\log$ \tabularnewline 
  Weibull (``\verb|weibull|'')                & $z_t>0$                  & $\verb|shape|>0$                 &  $\log$ \tabularnewline  
                                              &                          & $\verb|scale|>0$                 &  $\log$ \tabularnewline  
                                              &                          & $\verb|zero-mass|\in(0,1)$       &  $\text{logit}$ \tabularnewline 
  Wrapped Cauchy (``\verb|wrpcauchy|'')       & $z_t\in(-\pi,\pi]$       & $\verb|mean|\in(-\pi,\pi]$       &  $\tan(\verb|mean|/2)$ \tabularnewline  
                                              &                          & $\verb|concentration|\in(0,1)$   &  $\text{logit}$ \tabularnewline 
  \bottomrule
  \end{tabular}
  \footnotesize{$^1$Link functions $(g)$ relate natural scale parameters $({\boldsymbol \theta})$ to a $T \times K$ design matrix $({\mathbf X})$ and $K-$vector of working scale parameters $(\boldsymbol{\beta}\in \mathbb{R}^K)$ such that $g({\boldsymbol \theta})={\mathbf X}\boldsymbol{\beta}$. %In circular-circular regression models for the mean of angular distributions, we have $\theta_t = {\text atan2}( \sum_{k=1}^K z_{t,k} \sin(x_{t,k} \beta_k),1+\sum_{k=1}^K z_{t,k} \cos(x_{t,k} \beta_k))$, where $Z$ is a $T \times K$ matrix of positive real covariates (e.g. wind speed) and $X$ is a $T \times K$ matrix composed of angular covariates (e.g. wind direction).
}
\end{table}

%\subsection{Workflow}

\subsection{Data preparation and visualization}
For temporally-regular location data with negligible measurement error, the \verb|prepData| function is used to create a \verb|momentuHMMData| object that can be used for data visualization and further analysis. The arguments for \verb|prepData| include:
\begin{itemize}
  \item{\verb|data|} A data frame with $T+1$ rows including optionally a field `ID' (identifiers for different individuals), coordinates from which step length (`step') and turning angle (`angle') data streams are to be calculated, any additional data streams, and any covariates identified in the \verb|covNames| and \verb|angleCovs| arguments. Altervatively, \verb|data| can be a \verb|crwData| object returned by \verb|crawlWrap|.
  \item{\verb|type|} Coordinate type; `UTM' if easting-northing or `LL' if longitude-latitude.
  \item{\verb|coordNames|} Names of the two coordinate columns in \verb|data|. If \verb|coordNames=NULL| then step lengths, turning angles, and any location-based covariates (i.e., those specified by \verb|spatialCovs|, \verb|centers|, \verb|centroids|, and \verb|angleCovs|) are not calculated.
  \item{\verb|covNames|} Character vector indicating the names of any covariates in \verb|data|. Any variables in \verb|data| (other than ``ID'') that are not identified in \verb|covNames| or \verb|angleCovs| are assumed to be data streams.
  \item{\verb|spatialCovs|} List of Raster-class objects \citep{Hijmans2016} containing spatio-temporally referenced covariates. Covariates specified by \verb|spatialCovs| are extracted from the raster layer(s) based on the location data. Raster stacks may also be included, in which case the appropriate $z$ values (e.g. time, date) must also be included in \verb|data|.
  \item{\verb|centers|} 2-column matrix providing the coordinates for any activity centers (e.g., potential centers of attraction or repulsion) from which distance and angle covariates will be calculated based on the location data and returned in the \verb|momentuHMMData| object.
  \item{\verb|centroids|} List where each element is a data frame containing the x-coordinates ('x'), y-coordinates ('y'), and times for a centroid (i.e., a dynamic activity center for which the coordinates can change over time) from which distance and angle covariates will be calculated based on the location data and returned in the \verb|momentuHMMData| object.
  \item{\verb|angleCovs|} Character vector indicating the names of any circular-circular regression angular covariates in \verb|data| or \verb|spatialCovs| that need conversion from standard direction (in radians relative to the x-axis) to turning angle (relative to previous movement direction).
\end{itemize}
Summary plots of the \verb|momentuHMMData| object returned by \verb|prepData| can be created for any data stream or covariate using the generic \verb|plot| function.

If location data are temporally-irregular or subject to measurement error, then they are not suitable for \verb|prepData|. In this case, \verb|momentuHMM| can be used to perform a 2-stage multiple imputation approach \citep{McClintock2017}. We discuss this pragmatic approach to incorporating uncertainty attributable to observation error and temporal irreglarity into multivariate HMM analyses in section \ref{sec:mi}.

\subsection{HMM specification and fitting}
Once a \verb|momentuHMMData| object has been created using \verb|prepData|, then the data are ready to be passed to the generalized multivariate HMM-fitting function \verb|fitHMM|. There are many different options for specifying HMMs using \verb|fitHMM|, so here we will only focus on several of the most important and useful features (further details of all \verb|fitHMM| arguments are in the package documentation). The bare essentials of \verb|fitHMM| include the arguments:
\begin{itemize}
  \item{\verb|data|} A \verb|momentuHMMData| object
  \item{\verb|nbStates|} Number of latent states $(N)$
  \item{\verb|dist|} A named list indicating the probability distributions of the data streams.
  \item{\verb|estAngleMean|} An optional named list indicating whether or not to estimate the angle mean for data streams with angular distributions (e.g. turning angle). If not estimated (the default), the angle mean is fixed to 0.
  \item{\verb|formula|} Regression formula for the transition probability covariates
  \item{\verb|stationary|} Logical indicating whether or not the initial distribution is considered equal to the stationary distribution (must be \verb|FALSE| if \verb|formula| includes covariates) 
  \item{\verb|Par0|} A named list containing vectors of starting values for the state-dependent probability distribution parameters of each data stream
\end{itemize}
These seven arguments are all that are needed in order to fit the HMMs currently supported in \verb|moveHMM| \citep{MichelotEtAl2016}. For example, here is how the analysis of 15 ``wild haggis'' tracks described in \cite{MichelotEtAl2016} would be implemented using \verb|momentuHMM|:
<<load-library, echo=TRUE, eval=FALSE>>=
library(momentuHMM)
### Load raw data
rawHaggis<-read.csv("rawHaggises.csv")
### Process data
processedHaggis<-prepData(data=rawHaggis,covNames=c("slope","temp"))

### Fit HMM	
# initial step distribution natural scale parameters
stepPar0 <- c(1,5,0.5,3) # (mu_1,mu_2,sd_1,sd_2)
# initial angle distribution natural scale parameters 
anglePar0 <- c(0,0,1,8) # (mean_1,mean_2,concentration_1,concentration_2)       	
fitHaggis <- fitHMM(data = processedHaggis, nbStates = 2,
                    dist = list(step = "gamma", angle = "vm"),
                    Par0 = list(step = stepPar0, angle = anglePar0),
                    formula = ~ slope + I(slope^2),
                    estAngleMean = list(angle=TRUE))
@
Note that many of the arguments in \verb|fitHMM| are lists, with each element of the list corresponding to a data stream.  The list names provided in \verb|dist|, \verb|Par0|, and \verb|estAngleMean| (e.g. `step' and `angle') must therefore have a corresponding column in \verb|data| with the same name.  Additional data streams can be added to the model by simply adding the additional elements to these list arguments (see examples in sections \ref{sec:nfs} and \ref{sec:northernFulmar}).  

As seen above, the \verb|formula| argument can include many of the functions and operators commonly used to construct terms in R linear model formulas (e.g. \verb|a*b|, \verb|a:b|, \verb|cos(a)|). The \verb|formulaDelta| argument can be similarly used to specify covariate models for the initial distribution.  Unique to \verb|momentuHMM|, the \verb|formula| argument can also be used to specify transition probability matrix models that incorporate cyclical patterns (using the \verb|cosinor| special function; see example in section \ref{sec:elephant}), splines for explaining other more complicated patterns (e.g., \verb|bs| and \verb|ns| functions in the R base package \verb|splines|), and factor variables (e.g., \verb|formula=~ID| for individual-level effects).  By default the \verb|formula| argument applies to all state transition probabilities, but the special functions \verb|state|, \verb|toState|, and \verb|betaCol| allow for state- and parameter-specific formulas to be specified (see examples in sections \ref{sec:greySeal} and \ref{sec:northernFulmar}). While \verb|betaCol| allows a formula to be specified for a specific transition (e.g. state $3 \rightarrow 1$), \verb|state| and \verb|toState| allow a formula to be specified for all transitions from (e.g. $3 \rightarrow 1$, $3 \rightarrow 2$) and to (e.g. state $1 \rightarrow 3$, $2 \rightarrow 3$) specific states, respectively. The \verb|betaCons| argument allows for equality constraints among any of the transition probability parameters (e.g. $\gamma^{(t)}_{12}=\gamma^{(t)}_{21}$; see example in section \ref{sec:northernFulmar}). Specific state transition probabilities can also be fixed to zero (or any other value) using the \verb|fixPar| argument, which can be useful for incorporating more behavioral realism.  For example, \verb|fixPar| can be used to prohibit or enforce switching from one particular state to another (possibly as a function of spatio-temporal covariates). 

Similar to the \verb|formula| argument for state transition probability modelling, it is through the \verb|DM| argument of \verb|fitHMM| that models are specified for the state-dependent probability distribution parameters for each data stream.  \verb|DM| is a list argument containing an element for each data stream, but each element itself is also a list specifying the design matrix formulas for each parameter. For example, the following fits the exact same wild haggis model as above, but employs a user-specified (intercept-only) design matrix for the step length data stream:
<<fit-haggis, echo=TRUE, eval=FALSE>>=
stepDM <- list(mean = ~1, sd = ~1)

### Fit HMM	using user-specified DM
fitHaggisDM <- fitHMM(data = processedHaggis, nbStates = 2,
                      dist = list(step = "gamma", angle = "vm"),
                      DM = list(step = stepDM),
                      Par0 = list(step = log(stepPar0), angle = anglePar0),
                      formula = ~ slope + I(slope^2),
                      estAngleMean = list(angle=TRUE))
@
Note that when \verb|DM| is specified for a data stream, the initial parameter values (\verb|Par0|) for that data stream now correspond to columns of the resulting design matrix and must be on the working scale instead of the natural scale.  In this case, because the log link is used for the natural parameters of the gamma distribution, \verb|Par0$step| was specified on the log scale.  The functions \verb|getPar|, \verb|getPar0|, \verb|checkPar0|, and \verb|getParDM| are designed to assist users in the specification of design matrices and corresponding initial values on the working scale for any given model (see package documentation for further details).  \verb|DM| formulas are just as flexible as the \verb|formula| argument and, in addition to common linear model formula functions and operators, can also include cyclical cosinor models (see section \ref{sec:elephant}), splines, factor variables, and state-specific probability distribution parameter formulas (see examples in sections \ref{sec:turtle} and \ref{sec:greySeal}). As with the state transition probabilities, working parameters for probability distributions can also be fixed to user-specified values using the \verb|fixPar| argument.

Specification of design matrices using \verb|DM| is not limited to formulas. Alternatively, ``pseudo-design'' matrices can be specified, using an R matrix with rows corresponding to the natural parameters and columns corresponding to the working parameters. The elements in the matrix may be numeric or character strings containing model formula terms (see examples in sections \ref{sec:greySeal}, \ref{sec:harbourSeal}, and \ref{sec:northernFulmar}). Using a pseudo-design matrix for step length, the following is yet another way to implement the exact same wild haggis model:
<<fit-haggis-2, echo=TRUE, eval=FALSE>>=
stepDMp <- matrix(c(1,0,0,0,
                   0,1,0,0,
                   0,0,1,0,
                   0,0,0,1),4,4,byrow=TRUE)
rownames(stepDMp) <- c("mean_1","mean_2","sd_1","sd_2")
colnames(stepDMp) <- c("mean_1:(Intercept)","mean_2:(Intercept)",
                      "sd_1:(Intercept)","sd_2:(Intercept)")

### Fit HMM	using user-specified DM
fitHaggisDMp <- fitHMM(data = processedHaggis, nbStates = 2,
                       dist = list(step = "gamma", angle = "vm"),
                       DM = list(step = stepDMp),
                       Par0 = list(step = log(stepPar0), angle = anglePar0),
                       formula = ~ slope + I(slope^2),
                       estAngleMean = list(angle=TRUE))
@
\noindent (note that column and rownames for pseudo-design matrices are not required but can be useful). Pseudo-design matrices allow for the sharing of common working parameters (such as intercept terms) among natural scale parameters, and this can be used to constrain natural scale parameters (e.g., $\mu_1 \le \mu_2$) when used in tandem with the \verb|workBounds| argument (see sections \ref{sec:nfs}, \ref{sec:harbourSeal}, and \ref{sec:northernFulmar}).  This is particularly useful for preventing state label switching when repeatedly fitting the same HMM using multiple imputation methods (see section \ref{sec:mi}).

\subsection{Circular-circular regression model for the angle mean}
\label{sec:circ}
Another noteworthy \verb|fitHMM| argument, \verb|circularAngleMean|, is a list argument that enables users to specify circular-circular regression models for the mean $(\mu)$ parameter of angular distributions, such as the wrapped Cauchy and von Mises, instead of circular-linear models based on the tangent link function (Table \ref{tab:pdfs}). When \verb|circularAngleMean| is specified as \verb|TRUE| for any given angular data stream (e.g. turning angle), then a special link function based on \cite{RivestEtAl2016} is used:
\begin{equation}
  {\boldsymbol \mu}=\text{atan2}(\sin({\mathbf X}_\mu){\boldsymbol \beta}_\mu,1+\cos({\mathbf X}_\mu){\boldsymbol \beta}_\mu),
  \label{eq:circ}
\end{equation}
where ${\mathbf X}_\mu$ is a $T \times K$ matrix composed of the turning angles between $K$ angular covariates (e.g., wind direction, sea surface current direction) and the bearing of movement during the previous time step; that is, each element 
\begin{equation}
x_{t,k}=\text{atan2}(\sin(r_{t,k}-b_{t-1}),\cos(r_{t,k}-b_{t-1})) 
  \label{eq:angleCov}
\end{equation}
for angular covariate $r_{t,k}$ and $k=1,\ldots,K$ (note that \verb|prepData| and \verb|MIfitHMM| calculate ${\mathbf X}_\mu$ based on the \verb|angleCovs|, \verb|centers|, or \verb|centroids| arguments so users need not bother). Because this link function is designed for turning angles, a turning angle of 0 is provided as the reference angle (hence the ``$1+$'' preceeding the cosine term in Eq. \ref{eq:circ}).  Thus as a trade-off between biased and correlated movements, the working parameters $({\boldsymbol \beta}_\mu)$ for the expected turning angle at time $t$ weight the attractive (or repulsive) strengths of the angular covariates relative to directional persistence.  When all ${\boldsymbol \beta}_\mu=0$, the model reduces to a correlated random walk, but an increasingly biased random walk results as ${\boldsymbol \beta}_\mu$ gets larger (or smaller). Many interesting hypotheses about animal movmement can be addressed using circular-circular regression on movement direction, including the effects of wind, sea surface currents (see example in section \ref{sec:turtle}), centers of attraction or repulsion (see examples in sections \ref{sec:greySeal} and \ref{sec:northernFulmar}), group dynamic models (see example in section \ref{sec:groupModel}), and dynamic activity centers (see example in section \ref{sec:northernFulmar}.

The special function \verb|angleFormula| can be included in \verb|DM| formulas or pseudo-design matrices in order to model the circular-circular regression angle mean as a function of the relative strength (or importance) of angular covariates \citep{RivestEtAl2016}:
\begin{equation}
  {\boldsymbol \mu}=\text{atan2}(({\mathbf Z}_\mu \circ \sin({\mathbf X}_\mu)){\boldsymbol \beta}_\mu,1+({\mathbf Z}_\mu \circ \cos({\mathbf X}_\mu)){\boldsymbol \beta}_\mu),
  \label{eq:circ2}
\end{equation}
where ${\mathbf Z}_\mu$ is a $T \times K$ matrix of positive real covariates (e.g. wind speed, sea surface current speed) and $\circ$ is the Hadamard (i.e. element-wise) product.  The special function \verb|angleFormula| can also be used to specify group- or individual-level effects on the circular-circular regression angle mean coefficients $({\boldsymbol \beta}_\mu)$.

Also based on \cite{RivestEtAl2016}, the von Mises consensus distribution is a special von Mises circular-circular regression model where the concentration parameter $(\rho)$ depends on the level of agreement among short-term directional persistence (i.e. moving forward) and the angular covariates:
\begin{equation}
  {\boldsymbol \rho}=\kappa \sqrt{\left[({\mathbf Z}_\mu \circ \sin({\mathbf X}_\mu)){\boldsymbol \beta}_\mu\right]^2+\left[1+({\mathbf Z}_\mu \circ \cos({\mathbf X}_\mu)){\boldsymbol \beta}_\mu\right]^2}.
  \label{eq:consensus}
\end{equation}
Note that the von Mises consensus distribution is parameterized in terms of $\mu$ and $\kappa$ (see Table \ref{tab:pdfs}), but \verb|momentuHMM| returns and plots real parameter estimates in terms of $\mu$ and $\rho$. When all ${\boldsymbol \beta}_\mu$ are non-negative, then the minimum and maximum values for $\rho$ are $\kappa \lvert 1-\min({\mathbf Z}_\mu {\boldsymbol \beta}_\mu) \rvert$ and $\kappa \left[1+\max({\mathbf Z}_\mu {\boldsymbol \beta}_\mu)\right]$, respectively.  In the consensus model, $\kappa$ can be interpreted as the concentration towards a turning angle of zero (i.e. moving forward) when the angular covariate components perfectly cancel out. %For example, when $z_1=z_2$, $x_1 = -x_2$, and $\beta_1=\beta_2$:
%\begin{equation*}
%\rho = \kappa \sqrt{\left[z_1\sin(x_1)\beta_1 + z_2\sin(x_2)\beta_2\right]^2+\left[1+z_1\cos(x_1)\beta_1 + z_2\cos(x_2)\beta_2\right]^2} = \kappa.
%\end{equation*}
See section \ref{sec:turtle} for example code using \verb|angleFormula| and the von Mises consensus (``vmConsensus'') distribution.

\subsection{Multiple imputation}
\label{sec:mi}
When location data are temporally-irregular or subject to measurement error, then they are not suitable for standard maximum-likelihood HMM analyses based on the forward algorithm (Eq. \ref{eq:HMMlike}). In this case, \verb|momentuHMM| can be used to perform the 2-stage multiple imputation approach of \cite{McClintock2017}. The basic concept is to first employ a single-state (i.e., $N=1$) movement model that is relatively easy to fit but can accommodate location measurement error and temporally-irregular or missing observations \citep[e.g.][]{JohnsonEtAl2008}. The second stage involves repeatedly fitting the desired HMM to $m$ temporally-regular realizations of the position process drawn from the model output of the first stage.  Data streams or covariates that are dependent on location (e.g., step length, turning angle, habitat type, snow depth, sea surface temperature) will of course vary among the $m$ realizations of the position process, and the pooled inferences across the HMM analyses therefore reflect location uncertainty.  

There are three primary functions (\verb|MIfitHMM|, \verb|MIpool|, and \verb|crawlWrap|) for performing multiple imputation HMM analyses in \verb|momentuHMM|, and all rely on parallel processing to speed up computations. \verb|crawlWrap| is a wrapper function for fitting the continuous-time correlated random walk (CTCRW) model of \cite{JohnsonEtAl2008} to one or more tracks (subject to location measurement error and/or temporal irregularity) and then predicting temporally-regular tracks of the user's choosing (e.g. 15 min, hourly, daily) based on the CTCRW model output.  \verb|crawlWrap| returns a \verb|crwData| object that can be used to draw $m$ realization of the position process within the \verb|MIfitHMM| function.  \verb|MIfitHMM| is essentially a wrapper function for \verb|fitHMM| that repeatedly fits the same user-specified HMM to $m$ imputed data sets and stores the output from each of the $m$ model fits. If a \verb|crwData| object is provided, then \verb|MIfitHMM| will first draw $m$ imputations based on the \verb|crwData| output and then fit the specified HMM to each imputed data set. If users wish to use a movement model other than the CTCRW to account for measurement error and temporal irregularity \citep[e.g.][]{CalabreseEtAl2016,GurarieEtAl2017}, or if other observation error processes (e.g. missing data) are to be accounted for in the imputation step, \verb|MIfitHMM| can also be used for analysis of a list of $m$ \verb|momentuHMMData| objects that were imputed by the user. Based on the $m$ model fits, the \verb|MIpool| function calculates pooled estimates, standard errors, and confidence intervals for the working scale parameters, natural scale parameters (based on transformations of the pooled working parameters and mean or user-specified values for any covariates), state sequences, state probabilities, and activity budgets (i.e. the proportion of the $T$ times step assigned to each state) using standard multiple imputation formulae \citep{RubinSchenker1986,McClintock2017}.  \verb|MIpool| can be called separately or within \verb|MIfitHMM| (using the \verb|poolEstimates| argument), and the function returns a \verb|miSum| object containing the pooled output across all imputatons. See sections \ref{sec:nfs}, \ref{sec:turtle}, \ref{sec:greySeal}, and \ref{sec:harbourSeal} for example HMM analyses that use multiple imputation to account for location measurement error and temporally irregularity.

\subsection{Model visualization and diagnostics}
The generic \verb|plot| functions for \verb|momentuHMM| models (\verb|plot.momentuHMM| and \verb|plot.miSum|) plot the data stream histograms along with their corresponding estimated probability distributions, the estimated natural parameters and state transition probabilities as a function of any covariates included in the model, and the tracks of all individuals (color-coded by the most likely state sequence). By default, the probability distributions are plotted based on the means of any covariate values, but user-specified covariate values for the plots can be provided using the \verb|covs| argument.  When the argument \verb|plotCI=TRUE|, then confidence intervals for the natural parameters and state transition probabilities are also plotted. Confidence intervals are calculated from the working parameter estimates based on the delta method and finite-difference approximations of the first derivative for the transformation using the \verb|numDeriv::grad| function \citep{GilbertVaradhan2016}.  For multiple imputation analyses (\verb|plot.miSum|), all plots are based on the pooled parameter estimates and the means of any covariates (if not provided by the \verb|covs| argument) across each imputation. Using the argument \verb|errorEllipse|, \verb|plot.miSum| will include estimated location error ellipses in the plots of individual tracks. The functions \verb|plotSat|, \verb|plotSpatialCov|, and \verb|plotStates| (Table \ref{tab:functions}) provide further methods for visualizing model results.

Diagnostic tools include the calculation and plotting of pseudo-residuals \citep{ZucchiniEtAl2016} using the \verb|pseudoRes| and \verb|plotPR| functions, respectively. For discrete distributions (e.g.\ Bernoulli, Poisson), a continuity adjustment is used for calculating pseudo-residuals. Akaike's Information Criterion can be calculated for one or more models using the \verb|AIC.momentuHMM| function.

\subsection{Simulation}
The function \verb|simData| can be used to simulate multivariate HMM data from scratch or based on the estimated parameters of existing \verb|momentuHMM| or \verb|miSum| models.  The \verb|simData| arguments are very similar to those used for model specification in \verb|fitHMM| (e.g., \verb|dist|, \verb|DM|) and data preparation in \verb|prepData| (e.g., \verb|spatialCovs|, \verb|centers|), but they include additonal arguments, \verb|lambda| and \verb|errorEllipse|, for simulating location data subject to temporal irregularity and measurement error, respectively. The \verb|spatialCovs| argument allows for rasters of spatio-temporal covariate values to be utilized in simulation models, while the \verb|centers| argument allows activity centers to be incorporated. Thus \verb|simData| can be used to simulate more ecologically-realistic tracks (potentially subject to observation error) that can be useful for study design, power analyses, and assessing model performance. Goodness-of-fit can also be investigated by drawing simulated data sets from a fitted model and comparing them to observed properties of the data \citep{MoralesEtAl2004}.
  
\section{Examples}
\label{sec:example}
We will now demonstrate some of the capabilities of \verb|momentuHMM| using real telemetry data. These examples are intended for demonstration purposes only, and we do not claim these example analyses represent improvements relative to previous or alternative analyses for these data sets. While only some of the key workflow elements are included here, complete R code and further details for these analyses are available in the ``vignettes'' source directory.

<<load-vignette_inputs, results='hide', echo=FALSE, eval=TRUE, message=FALSE, warning=FALSE>>=
library(momentuHMM)
#example_wd <- ("~/Dropbox/vignette examples/")
example_wd <- ("~/Documents/Dropbox/current projects/moveHMM extension/momentuHMM/vignette examples/")

#source("vignette_examples.R")
load("vignette_inputs.RData")
@
\subsection{African elephant}
\label{sec:elephant}
As our first example, we use an African elephant ({\it Loxodonta africana}) bull track described in \cite{WallEtAl2014} and publicly available from the movebank.org data repository. 

We can load the data from the URL (this requires an Internet connection):
<<load-data, echo=TRUE, eval=TRUE>>=
URL <- paste0("https://www.datarepository.movebank.org/bitstream/handle/",
              "10255/move.373/Elliptical%20Time-Density%20Model%20%28Wall%",
              "20et%20al.%202014%29%20African%20Elephant%20Dataset%20%",
              "28Source-Save%20the%20Elephants%29.csv")
rawData <- read.csv(url(URL))
@ 

The data set contains two tracks; for this analysis, we only consider the first one. In addition to hourly locations, the tag also collected external temperature data. We subset the data frame to keep only the relevant rows and columns:
<<subset-data, echo=TRUE>>=
# select and rename relevant columns
rawData <- rawData[,c(11,3,4,5,6)]
colnames(rawData) <- c("ID","time","lon","lat","temp")

# only keep first track
rawData <- subset(rawData,ID==unique(ID)[1])
@

The data now has the following columns:
<<head-data, echo=TRUE>>=
head(rawData)
@ 

Location measurement error is negligible for these terrestrial GPS data, although about 1\% of the hourly observations collected between 22 March 2008 and 30 September 2010 are missing. Instead of simply ignoring these missing data, we can employ \verb|crawlWrap| to predict the missing locations based on the CTCRW model of \cite{JohnsonEtAl2008} prior to conducting our HMM analysis.

To use \verb|crawlWrap|, we convert times from factors to POSIX, and project the observed locations to UTM coordinates:
<<project-data, echo=TRUE, eval=FALSE>>=
# convert times from factors to POSIX
rawData$time <- as.POSIXct(rawData$time,tz="GMT")

# project to UTM coordinates using package rgdal
library(rgdal)
llcoord <- SpatialPoints(rawData[,3:4], 
                         proj4string=CRS("+proj=longlat +datum=WGS84"))
utmcoord <- spTransform(llcoord,CRS("+proj=utm +zone=30 ellps=WGS84"))

# add UTM locations to data frame
rawData$x <- attr(utmcoord,"coords")[,1]
rawData$y <- attr(utmcoord,"coords")[,2]
@ 

Then, we call \verb|crawlWrap| to fit a CTCRW model and predict hourly locations:
<<crawl, echo=TRUE, eval=FALSE>>=
# fit crawl model
crwOut <- crawlWrap(obsData=rawData, timeStep="hour",
                    theta=c(6.855, -0.007), fixPar=c(NA,NA))
@

Here the desired time step is specified by the \verb|timeStep| argument, and \verb|theta| and \verb|fixPar| arguments are the same as for \verb|crawl::crwMLE| \citep{Johnson2017}. For the choice of initial parameters in \verb|crawlWrap|, we refer the reader to the documentation of the package crawl, in particular \verb|crawl::crwMLE| and \verb|crawl::crwPredict|. We now have a complete set of temporally-regular location data.

Autocorrelation function (ACF) estimates suggest there are 24-hour cycles in the step length data, and this presents an opportunity to demonstrate the use of the \verb|cosinor| function for incorporating cyclical behavior in model parameters using \verb|momentuHMM|. We create a \verb|momentuHMMData| object, and the 24-hour cosinor model covariate: 
<<prepData, echo=TRUE, eval=FALSE>>=
# create momentuHMMData object from crwData object
elephantData <- prepData(data=crwOut, covNames="temp")

# add cosinor covariate based on hour of day
elephantData$hour <- as.integer(strftime(elephantData$time, format = "%H", tz="GMT"))
@

As seen here, the function \verb|prepData| can also be used for pre-processing the best predicted track data from \verb|crawlWrap| output. The 24-hour cosinor covariate (``hour'') is simply a set of integers $(0,1,\ldots,23)$ indicating the hour of day for each observation. The ACF plot of the step lengths, shown in Figure \ref{fig:elephantResults1}, was obtained with:
<<acf-before, echo=TRUE, eval=FALSE>>=
acf(elephantData$step[!is.na(elephantData$step)],lag.max=300)
@ 

Our aim is to fit a 2-state HMM to the elephant track that includes temperature effects on the turning angle concentration parameters and cycling temperature effects (with a 24-hour periodicity) on the step length and state transition probability parameters. Complex models such as this can require many parameters, and it can be challenging to choose good starting parameter values for the optimization. Here, we take an incremental approach, starting from a simpler model with no covariates. In \verb|momentuHMM|, the function \verb|getPar0| extracts initial parameters from a fitted (nested) HMM, given arguments for the more complex model.

For the covariate-free 2-state model, six initial parameters need to be chosen: for each state, the mean and standard deviation of the gamma distribution of step lengths, and the concentration of the wrapped Cauchy distribution of turning angles. Looking at the histograms of the step lengths and the turning angles (e.g.\ output by \verb|plot(elephantData)|) is often useful to choose good starting parameter values.

<<fitm1, echo=TRUE, eval=FALSE>>=
# label states
stateNames <- c("encamped","exploratory")
# distributions for observation processes
dist = list(step = "gamma", angle = "wrpcauchy")

# initial parameters
Par0_m1 <- list(step=c(100,500,100,200),angle=c(0.3,0.7))

# fit model
m1 <- fitHMM(data = elephantData, nbStates = 2, dist = dist, Par0 = Par0_m1, 
             estAngleMean = list(angle=FALSE), stateNames = stateNames)
@ 

To ensure convergence, we could also use the argument \verb|retryFits| to specify the number of attempts to minimize the negative log-likelihood based on random perturbations of the parameter estimates at the current minimum.\\

We can build on complexity, by including the temperature and time of day as covariates in the state transition probabilities. We use the function \verb|getPar0| to extract the new starting parameter values.
<<fitm2, echo=TRUE, eval=FALSE>>=
# formula for transition probabilities
formula <- ~ temp * cosinor(hour, period = 24)

# initial parameters (obtained from nested model m1)
Par0_m2 <- getPar0(model=m1, formula=formula)

# fit model
m2 <- fitHMM(data = elephantData, nbStates = 2, dist = dist, Par0 = Par0_m2$Par, 
             beta0=Par0_m2$beta, stateNames = stateNames, formula=formula)
@ 

The special function \verb|cosinor(hour, period = 24)| internally creates the cosinor model covariates, $\cos(2\pi \times hour / period)$ and $\sin(2\pi \times hour / period)$, and includes both terms (plus interactions with ``temp'') in the fitted model.\\

Finally, we can fit the more complex model, including the effect of temperature and time of day on the parameters of the state-dependent distributions of steps and angles.

<<fitm3, echo=TRUE, eval=FALSE>>=
# formulas for parameters of state-dependent observation distributions
DM <- list(step = list(mean = ~ temp * cosinor(hour, period = 24),
                       sd = ~ temp * cosinor(hour, period = 24)),
           angle = list(concentration = ~ temp))

# initial parameters (obtained from nested model m2)
Par0_m3 <- getPar0(model=m2, formula=formula, DM=DM)

# fit model
m3 <- fitHMM(data = elephantData, nbStates = 2, dist = dist, Par0 = Par0_m3$Par, 
             beta0 = Par0_m3$beta, DM = DM, stateNames = stateNames,
             formula = formula)
@ 

The above model \verb|m3| identifed a state of slow undirected movement (``encamped''), and a state of faster and more directed movement (``exploratory'') (Figure \ref{fig:elephantResults1}). For a fitted model, the function \verb|viterbi| computes the most likely state sequence:
<<viterbi, echo=TRUE, eval=FALSE>>=
# decode most likely state sequence
states <- viterbi(m3)
# derive percentage of time spent in each state
table(states)/nrow(elephantData)
@

Here, about \Sexpr{round(as.numeric(activityBudgets[1]),2)*100}\% of the steps were attributed to the ``encamped'' state, and \Sexpr{round(as.numeric(activityBudgets[2]),2)*100}\% were attributed to the ``exploratory'' state. 

We can use \verb|AIC(m1,m2,m3)| to compare the three fitted models in terms of AIC; here, \verb|m3| is overwhelmingly supported by the AIC when compared to alternative models with fewer covariates.

The model can be visualized with the generic function \verb|plot|, which was used for the plots shown in Figure \ref{fig:elephantResults2}, and the decoded track in Figure \ref{fig:elephantResults1}.

<<plotm3, echo=TRUE, eval=FALSE>>=
plot(m3, plotCI = TRUE, covs = data.frame(hour=12))
@ 

Interestingly, this model suggests step lengths and directional persistence for the ``encamped'' state decreased as temperature increased, step lengths for both states tended to decrease in the late evening and early morning, and transition probabilities from the ``encamped'' to ``exploratory'' state decreased as temperature increased (Figure \ref{fig:elephantResults2}). 

Model fit can be assessed using the pseudo-residuals, with the functions \verb|pseudoRes| and \verb|plotPR|. The residual ACF plot shown in Figure \ref{fig:elephantResults1} was produced by:
<<pseudoRes, echo=TRUE, eval=FALSE>>=
# compute pseudo-residuals for the steps and the angles
pr <- pseudoRes(m3)

# plot the ACF of step pseudo-residuals
acf(pr$stepRes[!is.na(pr$stepRes)],lag.max = 300)
@ 

Autocorrelation function plots of the pseudo-residuals indicate this model explained much of the periodicity in step length, although there does still appear to be some room for improvement.

\begin{figure}[htbp]
  \includegraphics[width=0.49\textwidth]{elephant_plotSat.png}
  \includegraphics[width=0.49\textwidth]{plot_elephantResults017.pdf}
  \includegraphics[width=0.49\textwidth]{plot_elephantResults015.pdf}
  \includegraphics[width=0.49\textwidth]{plot_elephantResults016.pdf}
  \caption{Plot of the elephant track produced using the `plotSat' function (top-left panel), autocorrelation function (ACF) plot of the corresponding step length data (top-right panel), plot of the Viterbi-decoded state sequence for the 2-state (``encamped'' and ``exploratory'') model generated using the generic `plot' function (bottom-left panel), and the step length pseudo-residual ACF plot for this model using the `plotPR' function (bottom-right panel).}
  \label{fig:elephantResults1}
\end{figure}

\begin{figure}[htbp]
  \centering
  \includegraphics[width=0.35\textwidth]{plot_elephantResults009.pdf}
  \includegraphics[width=0.35\textwidth]{plot_elephantResults012.pdf} \\
  \includegraphics[width=0.35\textwidth]{plot_elephantResults001.pdf} 
  \includegraphics[width=0.35\textwidth]{plot_elephantResults002.pdf} \\
  \includegraphics[width=0.35\textwidth]{plot_elephantResults010.pdf}
  \includegraphics[width=0.35\textwidth]{plot_elephantResults013.pdf} \\
  \caption{Selected plots for the 2-state (``encamped'' and ``exploratory'') African elephant example generated using the generic 'plot' function. Top panels present histograms of the step length (top-left) and turning angle (top-right) data along with the estimated state-dependent probability distributions based on the mean temperature (temp = 29.7 degrees celsius) at 12:00 GMT (hour = 12). Middle panels present estimates (and 95\% confidence intervals) for the step length mean parameter of the ``encamped'' state as a function of temperature and hour of day.  Bottom-left panel presents estimates for the turning angle concentration parameter of the ``encamped'' state as a function of temperature.  Bottom-right panel presents estimated state transition probabilities (1 = ``encamped'', 2 = ``exploratory'') as a function of temperature at 12:00 GMT.}
  \label{fig:elephantResults2}
\end{figure}

\subsection{Northern fur seal}
\label{sec:nfs}
In our second example, we use the northern fur seal ({\it Callorhinus ursinus}) example from \cite{McClintockEtAl2014b} to demonstrate the use of additional data streams for distinguishing behaviors with similar horizontal trajectories in a multivariate HMM. The data consist of 241 temporally-irregular Fastloc GPS locations obtained during a foraging trip of a nursing female near the Pribilof Islands of Alaska, USA, from 10-17 October 2007. The tag included time-depth recording capabilities, and the dive activity data were summarized as the number of foraging dives over $T=228$ temporally-regular 1 hr time steps. To fit the $N=3$ state (1=``resting'', 2=``foraging'', 3=``transit'') of \cite{McClintockEtAl2014b} using \verb|momentuHMM|, we first used \verb|crawlWrap| to predict temporally-regular locations at 1 hr time steps assuming a bivariate normal measurement error model and merged the results with the foraging dive data using the \verb|crawlMerge| function. Then multiple imputation was used to account for location measurement error by repeatedly fitting the HMM to \verb|nSims| realizations of the position process using \verb|MIfitHMM|:% in parallel on \verb|ncores| cores of a (multi-core) processor:
<<fit-nfs, echo=TRUE, eval=FALSE>>=
nbStates <- 3
stateNames <- c("resting", "foraging", "transit")
dist <- list(step = "gamma", angle = "wrpcauchy", dive = "pois")
Par0 <- getParDM(nbStates = nbStates, dist = dist,
                 Par = Par, DM = DM, workBounds = workBounds,
                 estAngleMean = list(angle = FALSE))
Fixpar <- list(dive = c(-100, NA, NA))
nfsFits <- MIfitHMM(crwOut, nSims = 100, nbStates = nbStates, dist = dist,
                    Par0 = Par0, DM = DM, workBounds = workBounds,
                    estAngleMean = list(angle = FALSE), 
                    fixPar = fixPar, retryFits = 30,
                    stateNames=stateNames)
plot(nfsFits)
@
Here we specified a gamma distribution for step length (`step'), wrapped Cauchy distribution for turning angle (`angle'), and Poisson distribution for the number of foraging dives (`dive'). The function \verb|getParDM| was used to organize the starting values for the data stream working parameters (\verb|Par0|) in the correct format based on \verb|DM|, \verb|workBounds|, and estimates of the natural parameters (\verb|Par|) from \cite{McClintockEtAl2014b}. The \verb|DM| and \verb|workBounds| arguments were specified to avoid label switching among the \verb|nSims| imputed data model fits and enforce similar state-dependent probability distribution constraints as \cite{McClintockEtAl2014b}; for example, constraining the Poisson rate parameters such that the ``foraging'' state tends to have higher numbers of foraging dives than the ``transit'' state ($\lambda_2 > \lambda_3$; see Eq. \ref{eq:link} in section \ref{sec:harbourSeal} for more details on parameter constraints using \verb|DM| in conjunction with the \verb|userBounds| and \verb|workBounds| arguments). To prohibit foraging dives for the ``resting'' state, we used the \verb|fixPar| argument to effectively fix the Poisson rate parameter to zero on the natural scale (i.e. $\lambda_1 \approx 0$).  To help deal with the problem of convergence to local maxima, the \verb|retryFits| argument allows users to specify the number of times to attempt to re-fit each model using random perturbations of the parameter estimates as the starting values for optimization.

The results are very similar to those of the discrete-time model of \cite{McClintockEtAl2014b}, with periods of foraging often followed by resting (Figure \ref{fig:nfsResults}).  The ``activity budgets'' (i.e. the proportion of time steps allocated to each state) calculated by \verb|MIpool| based on the estimated state sequences for each imputation were \Sexpr{round(as.numeric(nfsTimeInStates$est["resting"]),2)} (95\% CI: \Sexpr{round(as.numeric(nfsTimeInStates$lower["resting"]),2)}$-$\Sexpr{round(as.numeric(nfsTimeInStates$upper["resting"]),2)}) for ``resting'', \Sexpr{round(as.numeric(nfsTimeInStates$est["foraging"]),2)} (95\% CI: \Sexpr{round(as.numeric(nfsTimeInStates$lower["foraging"]),2)}$-$\Sexpr{round(as.numeric(nfsTimeInStates$upper["foraging"]),2)}) for ``foraging'', and \Sexpr{round(as.numeric(nfsTimeInStates$est["transit"]),2)} (95\% CI: \Sexpr{round(as.numeric(nfsTimeInStates$lower["transit"]),2)}$-$\Sexpr{round(as.numeric(nfsTimeInStates$upper["transit"]),2)}) for ``transit''.

\begin{figure}[htbp]
  \includegraphics[width=0.49\textwidth,page=1]{plot_nfsResults}
  \includegraphics[width=0.49\textwidth,page=2]{plot_nfsResults}
  \includegraphics[width=0.49\textwidth,page=3]{plot_nfsResults}
  \includegraphics[width=0.49\textwidth,page=4]{plot_nfsResults}
  \caption{Plots of the northern fur seal example results generated using the generic `plot' function. The estimated probability distributions for step length (top-left panel), turning angle (top-right panel), and number of foraging dives (bottom-left panel) for the 3-state (``resting'', ``foraging'', and ``transit'') model are plotted along with histograms of these data streams. The temporally-regular predicted locations (and 95\% ellipsoidal confidence bands) and estimated states are plotted in the bottom-right panel. All estimates are pooled across multiple imputations of the position process and thus reflect uncertainty attributable to location measurement error and temporally-irregular observations.}
  \label{fig:nfsResults}
\end{figure}
\subsection{Loggerhead turtle}
\label{sec:turtle}
For our third example, we demonstrate how to model movement direction and step length as a function of angular covariates using hitherto unpublished loggerhead turtle ({\it Caretta caretta}) data for a captive-raised juvenile released in 2012 on the coast of North Carolina, USA. The data consist of 165 temporally-irregular Argos locations subject to measurement error and rasters of daily ocean surface currents collected between 20 November and 19 December 2012. Assuming a gamma distribution for step length $(l_t)$ and a wrapped Cauchy distribution for turning angle $(\phi_t)$, we model the mean step length parameter $(\mu^l_t)$ as a function of ocean surface current speed $(w_t)$ and direction $(r_t)$ relative to the bearing of movement $(b_t)$:
\begin{equation}
  \mu^l_t=\exp(\beta^l_0+\beta^l_1 w_t \cos(b_t-r_t)),
  \label{eq:turtleMeanStep}
\end{equation}
%or, equivalently,
%\begin{equation}
%  \mu^l_t=\exp(\beta^l_0+\beta^l_1 w_t \cos(\phi_t-d_t))
%  \label{eq:turtleMeanStep}
%\end{equation}
and the turning angle mean parameter $(\mu^\phi_t)$ as a trade-off between short-term directional persistence and bias in the direction of ocean surface currents using the circular-circular regression link function:
\begin{equation}
  \mu^\phi_t=\text{atan2}(\sin(d_t) \beta^\phi,1+\cos(d_t)\beta^\phi),
    \label{eq:turtleMeanAngle}
\end{equation}
where $d_t=\text{atan2}(\sin(r_t-b_{t-1}),\cos(r_t-b_{t-1}))$.

We wish to fit a 2-state HMM to the turtle data, with a ``foraging'' state unaffected by currents and a ``transit'' state potentially influenced by ocean surface currents as in Eqs. \ref{eq:turtleMeanStep} and \ref{eq:turtleMeanAngle}.  We used \verb|crawlWrap| to predict $T=350$ temporally-regular locations at 2 hr time steps assuming a bivariate normal measurement error model that accounts for the Argos location quality class  (i.e. 3,2,1,0,A,B) of each observation. We then again used multiple imputation to account for location uncertainty by repeatedly fitting the HMM to \verb|nSims| realizations of the position process using \verb|MIfitHMM|. We first draw \verb|nSims| realizations of the position process and extract the corresponding spatial covariates from the raster bricks for ocean surface current speed (``speedBrick'') and direction (``dirBrick'') using \verb|MIfitHMM| with \verb|fit=FALSE|:
<<fit-turtle, echo=TRUE, eval=FALSE>>=
miTurtleData <- MIfitHMM(crwOut, nSims = 100, fit=FALSE,
                 spatialCovs = list(w = speedBrick, d = dirBrick, r = dirBrick),
                 angleCovs = "d")
@
When the \verb|fit| argument is \verb|FALSE|, \verb|MIfitHMM| returns a list of length \verb|nSims| composed of \verb|momentuHMMData| objects (\verb|miData|). For convenience and ease of interpretation, we manually added an additional covariate $\left({\text angle\_osc}=\cos(b_t-r_t)\right)$ to each of the imputed data sets and fitted the 2-state HMM using Eqs. \ref{eq:turtleMeanStep} and \ref{eq:turtleMeanAngle} for state 2 (``transit''):
%<<prep-turtle, echo=TRUE, eval=FALSE>>=
%for(j in 1:30){
%  miTurtleData$miData[[j]]$angle_osc <- miTurtleData$miData[[j]]$angle
%                                        -miTurtleData$miData[[j]]$osc_dir
%}
%@
%Now the imputed data are ready to be fitted to the 2-state HMM using Eqs. \ref{eq:turtleMeanStep} and \ref{eq:turtleMeanAngle} for state 2 (``transit''):
<<fit-turtle-2, echo=TRUE, eval=FALSE>>=
nbStates<-2
dist <- list(step = "gamma", angle = "wrpcauchy")
DM <- list(step = list(mean = ~state2(w:angle_osc), sd = ~1),
           angle = list(mean = ~state2(d), concentration= ~1))
turtleFits <- MIfitHMM(miTurtleData$miData, nbStates = nbStates, dist = dist, 
                       Par0 = Par0, DM = DM, 
                       estAngleMean = list(angle = TRUE),
                       circularAngleMean = list(angle = TRUE))
plot(turtleFits, plotCI = TRUE, covs = data.frame(angle_osc = cos(0)))
@
Note that the \verb|state2| special function in \verb|DM| indicates the covariate formulas are specific to state 2 (``transit'') and the \verb|circularAngleMean| argument indicates that circular-circular regression link function is to be used on the mean turning angle parameter as in Eq. \ref{eq:turtleMeanAngle}. 

For the ``transit'' state, pooled parameter estimates indicated step lengths increased with ocean surface current speed and as the bearing of movement aligned with ocean surface current direction ($\beta^l_1=\Sexpr{round(as.numeric(turtle_miSum$Par$beta$step$est[,"mean_2:w:angle_osc"]),2)}, \text{95\% CI: } \Sexpr{round(as.numeric(turtle_miSum$Par$beta$step$lower[,"mean_2:w:angle_osc"]),2)}-\Sexpr{round(as.numeric(turtle_miSum$Par$beta$step$upper[,"mean_2:w:angle_osc"]),2)}$; Figure \ref{fig:turtleResults}). The estimated wrapped Cauchy distribution for turning angle had mean angles $(\mu^\phi_t)$ biased towards the direction of ocean surface currents for each time step $(\beta^\phi=\Sexpr{round(as.numeric(turtle_miSum$Par$beta$angle$est[,"mean_2:(d)"]),2)}, \text{95\% CI: } \Sexpr{round(as.numeric(turtle_miSum$Par$beta$angle$lower[,"mean_2:(d)"]),2)}-\Sexpr{round(as.numeric(turtle_miSum$Par$beta$angle$upper[,"mean_2:(d)"]),2)})$, with concentration parameter $\rho^\phi_2=\Sexpr{round(as.numeric(turtle_miSum$Par$real$angle$est["concentration","state 2"]),2)}$ (95\% CI: \Sexpr{round(as.numeric(turtle_miSum$Par$real$angle$lower["concentration","state 2"]),2)}$-$\Sexpr{round(as.numeric(turtle_miSum$Par$real$angle$upper["concentration","state 2"]),2)}) indicating turning angles were concentrated at $\mu^\phi_t$. Thus movement during the ``transit'' state appears to strongly follow ocean surface currents (mean ${\text angle\_osc}=\Sexpr{round(mean(turtle_miSum$data$angle_osc[which(turtle_miSum$Par$states==2)]),2)},{\text SD}=\Sexpr{round(sd(turtle_miSum$data$angle_osc[which(turtle_miSum$Par$states==2)]),2)}$), while movement during the ``foraging'' state exhibited shorter step lengths $(\mu^l_1=\Sexpr{round(as.numeric(turtle_miSum$Par$real$step$est["mean","state 1"]),0)}{\text m}, \text{95\% CI: } \Sexpr{round(as.numeric(turtle_miSum$Par$real$step$lower["mean","state 1"]),0)}-\Sexpr{round(as.numeric(turtle_miSum$Par$real$step$upper["mean","state 1"]),0)})$ perpendicular to ocean surface currents (mean ${\text angle\_osc}=\Sexpr{round(mean(turtle_miSum$data$angle_osc[which(turtle_miSum$Par$states==1)]),2)},{\text SD}=\Sexpr{round(sd(turtle_miSum$data$angle_osc[which(turtle_miSum$Par$states==1)]),2)}$), with some directional persistence $(\rho^\phi_1=\Sexpr{round(as.numeric(turtle_miSum$Par$real$angle$est["concentration","state 1"]),2)}, \text{95\% CI: } \Sexpr{round(as.numeric(turtle_miSum$Par$real$angle$lower["concentration","state 1"]),2)}-\Sexpr{round(as.numeric(turtle_miSum$Par$real$angle$upper["concentration","state 1"]),2)})$. The turtle spent \Sexpr{round(as.numeric(turtle_miSum$Par$timeInStates$est["state 1"]),2)} (95\% CI: \Sexpr{round(as.numeric(turtle_miSum$Par$timeInStates$lower["state 1"]),2)}$-$\Sexpr{round(as.numeric(turtle_miSum$Par$timeInStates$upper["state 1"]),2)}) of the 2 hr time steps in the ``foraging'' state and \Sexpr{round(as.numeric(turtle_miSum$Par$timeInStates$est["state 2"]),2)} (95\% CI: \Sexpr{round(as.numeric(turtle_miSum$Par$timeInStates$lower["state 2"]),2)}$-$\Sexpr{round(as.numeric(turtle_miSum$Par$timeInStates$upper["state 2"]),2)}) of time steps in the ``transit'' state as it travelled northeast along a predominant current until it (presumably) found an attractive foraging patch (Figure \ref{fig:turtleResults}).

\begin{figure}[htbp]
  \centering
  \includegraphics[width=0.32\textwidth]{plot_turtleResults001.pdf}
  \includegraphics[width=0.32\textwidth]{plot_turtleResults002.pdf}
  \includegraphics[width=0.32\textwidth]{plot_turtleResults004.pdf}
  %\begin{adjustbox}{trim=0cm 0.25cm 0cm 1.5cm}
    \includegraphics[width=\textwidth]{plot_turtleResults2.pdf}
  %\end{adjustbox}
  \caption{Selected results from the loggerhead turtle example. Top panels include estimates and 95\% confidence intervals for the mean step length parameter as a function of ocean surface current speed $(w)$ when ocean surface current direction $(r_t)$ is the same as the bearing $(b_t)$ of movement (i.e. ${\text angle\_osc}=\cos(b_t-r_t)=1$; top-left panel), mean step length parameter as a function of ${\text angle\_osc}$ at the mean ocean surface current speed ($w=0.46$ m/s; top-middle panel), and mean turning angle parameter as a function of $d_t=\text{atan2}(\sin(r_t-b_{t-1}),\cos(r_t-b_{t-1}))$ (top-right panel). Bottom panel plots the pooled track, 95\% error ellipse confidence bands, and state (orange = ``foraging'', blue = ``transit'') estimates based on multiple imputations of the position process relative to ocean surface current speed (m/s) and direction on 2 December 2012.}
  \label{fig:turtleResults}
\end{figure}

It may often make more sense to weight angular covariates (such as ocean surface current direction) by their relative strength or importance. For example, weak ocean surface currents may be less likely to influence movement direction than strong ocean surface currents. This could easily be included in our turtle model using the \verb|angleFormula(cov, strength, by)| special function in \verb|DM|, where \verb|cov| is an angle covariate (e.g. wind direction), \verb|strength| is an optional positive real covariate (e.g. wind speed), and \verb|by| is an optional factor variable for individual- or group-level effects (e.g. ID, sex):
<<fit-turtle-3, echo=TRUE, eval=FALSE>>=
DM$angle = list(mean = ~state2(angleFormula(d, strength = w)), 
                concentration= ~1))
@
\noindent which would yield the following model for the ``transit'' state mean angle parameter:
\begin{equation}
  \mu^\phi_t=\text{atan2}(w_t \sin(d_t) \beta^\phi,1+w_t \cos(d_t)\beta^\phi).
    \label{eq:turtleMeanAngle2}
\end{equation}

Still another option would be to use the von Mises consensus model, where the concentration parameter would now depend on the level of agreement between short-term directional persistence (i.e. going forward) and ocean surface currents: 
<<fit-turtle-4, echo=TRUE, eval=FALSE>>=
dist$angle = "vmConsensus"
DM$angle = list(mean = ~state2(angleFormula(d, strength = w)), 
                kappa = ~1)
@
\noindent which would yield the following model for the ``transit'' state concentration parameter:
\begin{equation}
  \rho^\phi_t = \kappa \sqrt{[w_t \sin(d_t) \beta^\phi]^2+[1+w_t \cos(d_t)\beta^\phi]^2}.
    \label{eq:turtleMeanAngle3}
\end{equation}
If there were multiple turtles in this dataset, then individual-level effects could be included on $\mu^\phi$ by simply specifying \verb|angleFormula(d, strength = w, by = ID)| or \verb|angleFormula(d, by = ID)| (with no strength effects).


\subsection{Grey seal}
\label{sec:greySeal}
For our next example, we perform a similar analysis of a grey seal ({\it Halichoerus grypus}) track that was originally conducted by \cite{McClintockEtAl2012} using Bayesian methods and (computationally-intensive) Markov chain Monte Carlo. The data consist of 1045 temporally-irregular Fastloc GPS locations collected in the North Sea between 9 April and 11 August 2008. Because the seal repeatedly visited the same haul-out and foraging locations, it provides a nice example for demonstrating how to implement biased movements relative to activity centers using \verb|momentuHMM|. \cite{McClintockEtAl2012} fitted a 5-state model to these data that included three center of attraction states, with movement biased towards two haul-out sites (``Abertay'' and ``Farne Islands'') and a foraging area (``Dogger Bank''), and two ``exploratory'' states (''low speed'', ''high speed'') that were unassociated with an activity center. After using \verb|crawlWrap| to predict $T=1515$ temporally-regular locations at 2 hr time steps including a bivariate normal measurement error model, we can perform a very similar analysis to \cite{McClintockEtAl2012} in \verb|momentuHMM| by using the \verb|centers| argument and state-specific functions for the probability distribution parameters. A cluster analysis on the observed locations using the R package \verb|dtwclust| \citep{SardaEspinosa2017} identified three centroids with coordinates that were nearly identical to the three activity centers (``Abertay'', ``Farne Islands'', and ``Dogger Bank'') identified by \cite{McClintockEtAl2012}. We use these coordinates to derive covariates relative to the activity centers when drawing \verb|nSims| realizations of the position process:
<<fit-grey-seal, echo=TRUE, eval=FALSE>>=
crwSim <- MIfitHMM(crwOut, nSims = 100, fit=FALSE,
                  center = centers)
@
Specifying the \verb|centers| argument results in the calculation of two covariates for each activity center: the distance (with `.dist' suffix) and angle (with `.angle' suffix) from each location at time $t$. These covariates can then be used to model parameters as a function of the distance and angle to activity centers for each time step:
<<spec-grey-seal, echo=TRUE, eval=FALSE>>=
dist <- list(step = "weibull", angle = "wrpcauchy")
distFormula <- ~state1(I(Abertay.dist>2500)) + state2(I(Farne.dist>2500)) 
                  + state3(I(Dogger.dist>15000))
angleFormula <- ~state1(Abertay.angle) + state2(Farne.angle) 
                  + state3(Dogger.angle)
stepDM <- list(shape = distFormula, scale = distFormula)
angleDM <- list(mean = angleFormula, concentration = distFormula)
DM <- list(step = stepDM, angle = angleDM)
@
Similar to \cite{McClintockEtAl2012}, we assume a Weibull distribution for step length where both the shape and scale parameter depend on the distance from the location at time $t$ to each activity center. For the activity centers on land (``Abertay'' and ``Farne''), we allow the (state-dependent) step length parameters to change when the seal is beyond 2500m of the haulout. For the ``Dogger'' activity center, we allow the parameters to change when the seal is beyond 15000m of this (presumably) foraging area. We thus allow the movement behavior to change within these activity center states upon entering or leaving the vicinity of these sites.  We assume a wrapped Cauchy distribution for turning angle with (state-dependent) mean angle derived from the direction to each activity center at time $t$, and the concentration parameter is modeled similarly to the step length parameters. For the two ``exploratory'' states, we assumed they are simple random walks unaffected by proximity to activity centers. To complete our model specification, we use the \verb|knownStates| argument to assign the seal to the corresponding activity center state whenever it was within the 2500m (haul-out area) or 15000m (foraging area) thresholds for each imputed data set:
<<fit-grey-seal-2, echo=TRUE, eval=FALSE>>=
greySealFits <- MIfitHMM(miDat, nSims = 400,
                         nbStates = 5, dist = dist,
                         Par0 = Par0, beta0 = beta0, fixPar = fixPar,
                         formula = distFormula,
                         estAngleMean = list(angle=TRUE), 
                         circularAngleMean = list(angle=TRUE),
                         DM = DM, knownStates = knownStates)
plot(greySealFits, plotCI = TRUE)
@
As with the step length and turning angle concentration parameters, the state transition probabilities are also allowed to change as a function of distance to activity centers (as specified by the \verb|formula| argument). The starting values (\verb|Par0| and \verb|beta0|) for each imputation were extracted from a single HMM fitted to the best predicted locations from \verb|crawlWrap|, and \verb|fixPar| was used to remove short-term directional persistence (and thus formulate the model as a mixture of biased and simple random walks). 

Estimated activity budgets for the 5 states of this multiple imputation HMM were \Sexpr{round(greySealTimeInStates$est["Abertay"],2)} $(\Sexpr{round(greySealTimeInStates$lower["Abertay"],2)}-\Sexpr{round(greySealTimeInStates$upper["Abertay"],2)})$ for the ``Abertay'' haul-out state, \Sexpr{round(greySealTimeInStates$est["Farne"],2)} $(\Sexpr{round(greySealTimeInStates$lower["Farne"],2)}-\Sexpr{round(greySealTimeInStates$upper["Farne"],2)})$ for the ``Farne Islands'' haul-out state, \Sexpr{round(greySealTimeInStates$est["Dogger"],2)} $(\Sexpr{round(greySealTimeInStates$lower["Dogger"],2)}-\Sexpr{round(greySealTimeInStates$upper["Dogger"],2)})$ for the ``Dogger Bank'' foraging state, \Sexpr{round(greySealTimeInStates$est["low"],2)} $(\Sexpr{round(greySealTimeInStates$lower["low"],2)}-\Sexpr{round(greySealTimeInStates$upper["low"],2)})$ for a low-speed ``exploratory'' state, and \Sexpr{round(greySealTimeInStates$est["high"],2)} $(\Sexpr{round(greySealTimeInStates$lower["high"],2)}-\Sexpr{round(greySealTimeInStates$upper["high"],2)})$ for a high-speed ``exploratory'' state. All three activity center states exhibited shorter step lengths and less biased movements when within the vicinity of these targets (Figure \ref{fig:greySealResults}). Results from this analysis were thus very similar to those of \cite{McClintockEtAl2012}, but this implementation required far less computation time and no custom model-fitting algorithms. 

The \verb|simData| function can be used to simulate tracks from a fitted model:
<<sim-grey-seal, echo=TRUE, eval=FALSE>>=
greySealSim<-simData(model = greySealFits, centers = centers,
                     initialPosition = centers[1,],
                     obsPerAnimal = 1515)
@
A simulated track is presented along with the fitted track in Figure \ref{fig:greySealStateSims}. While potentially useful for study design, power analysis, and prediction, the \verb|simData| function can also be helpful in assessing goodness of fit by repeatedly drawing simulated data sets from a fitted model and comparing them to observed properties of the data \citep[e.g.][]{MoralesEtAl2004}.
\begin{figure}[htbp]
  \includegraphics[width=0.49\textwidth]{plot_greySealResults002.pdf}
  \includegraphics[width=0.49\textwidth]{plot_greySealResults009.pdf}\\
  \includegraphics[width=0.49\textwidth]{plot_greySealResults006.pdf}
  \includegraphics[width=0.49\textwidth]{plot_greySealResults013.pdf}\\
  \caption{Selected results from the grey seal example. Panels include estimates and 95\% confidence intervals for the ``Abertay'' haul-out state step length scale parameter as a function of distance in meters (`Abertay.dist'; top-left panel), ``Abertay'' haul-out state turning angle concentration parameter as a function of distance (top-right panel), ``Dogger Bank'' foraging state step length scale parameter as a function of distance (`Dogger.dist'; bottom-left panel), and the ``Dogger Bank'' foraging state turning angle concentration parameter as a function of distance (bottom-right panel).}
  \label{fig:greySealResults}
\end{figure}

\begin{figure}[htbp]
  \centering
  %\begin{adjustbox}{trim=0cm 0.25cm 0cm 1.5cm}
    \includegraphics[width=0.8\textwidth]{plot_greySealResults1.png}\\
    \includegraphics[width=0.8\textwidth]{plot_greySealResults2.png}
  %\end{adjustbox}
  \caption{Fitted and simulated tracks from the grey seal example. This seal tended to move in a clockwise fashion between two haul-out locations (``Abertay'' and ``Farne Islands'') and a foraging area (``Dogger Bank'') in the North Sea. Top panel plots the pooled track, 95\% error ellipse confidence bands, and state estimates based on the 5-state HMM fitted to multiple imputations of the position process. Red points indicate the locations of the three activity centers. Black points indicate the (temporally-irregular) observed locations. Bottom panel presents the locations and states for a track simulated from the fitted model using the `simData' function.}
  \label{fig:greySealStateSims}
\end{figure}

\subsection{Southern elephant seal}
Here, we analyse the southern elephant seal (\emph{Mirounga leonina}) data from \cite{MichelotEtAl2017} using \verb|momentuHMM|. The data set consists of 15 tracks, each encompassing (at most) one foraging trip, starting from Kerguelen Island. We want to fit the model described by \cite{MichelotEtAl2017}, with the four following states:
\begin{enumerate}
\item outbound trip from the colony to the ice;
\item searching;
\item foraging;
\item inbound trip from the ice to the colony.
\end{enumerate}

The data set has three columns: ``ID'' (track ID), ``x'' (longitude), and ``y'' (latitude):

<<head-ses, echo=TRUE, eval=TRUE>>=
head(tracks)
@ 

From the locations, we use \verb|prepData| to derive the step lengths and turning angles, as well as the distance and bearing (relative to previous movement direction as in Eq. \ref{eq:angleCov}) to the Kerguelen Island colony (with coordinates 70$^{\circ}$ longitude and -49$^{\circ}$ latitude):
<<prep-ses, echo=TRUE>>=
center <- matrix(c(70,-49),nrow=1,dimnames=list("colony"))
data <- prepData(data=tracks, type="LL", centers=center)
@
Note that distances are in kilometers and angles are based on initial bearings \citep[using geosphere::bearing;][]{Hijmans2016geo} when calculated from longitude and latitude coordinates.  

\subsubsection{Model 1: no covariates}
We start by fitting a covariate-free 4-state correlated random walk model, which we will use to extract starting parameter values for more complex models. We use the argument \verb|fixPar| to fix some transition probabilities to zero, following \cite{MichelotEtAl2017}.  We set to \verb|NA| the columns of unconstrained transition probabilities, and we fix the intercept of the other columns to a large negative number (here $-100$) to set the corresponding transition probabilities to be virtually zero (i.e.\ impossible transition).  As in \cite{MichelotEtAl2017}, we set transition probabilities from outbound to forage, outbound to inbound, search to outbound, forage to outbound, forage to inbound, inbound to outbound, inbound to search, and inbound to forage to be effectively zero.
<<m1-ses, echo=TRUE, eval=FALSE>>=
stateNames <- c("outbound","search","forage","inbound")

# initial parameters
stepPar0 <- c(25,5,1,25,10,5,3,10)
anglePar0 <- c(15,5,2,15)

# constrain transition probabilities
fixbeta <- matrix(c(NA,-100,-100,-100,NA,NA,-100,NA,-100,-100,-100,-100),
                  nrow=1)

m1 <- fitHMM(data=data, nbStates=4, dist=list(step="gamma",angle="vm"), 
             Par0=list(step=stepPar0, angle=anglePar0),
             fixPar=list(beta=fixbeta), stateNames = stateNames)
@ 

\subsubsection{Model 2}
This model mimics the formulation of \cite{MichelotEtAl2017}. We model the effect of the distance to colony on the transition probability from outbound to search, and of the time since departure on the transition probability from search to inbound.
<<formula-ses, echo=TRUE, eval=FALSE>>=
# time spent since left colony
time <- NULL
for(id in unique(data$ID)) {
    nbSubObs <- length(which(data$ID==id))
    
    # approximately in months for interval = 9.6h
    time <- c(time, (1:nbSubObs)/75)
}

data$time <- time

# compute time since departure and include in formula below
formula <- ~ colony.dist + time
@

As before, we constrain the transition probability matrix to prevent some of the transitions (e.g.\ from forage to inbound, etc.). We define a $3 \times 12$ matrix for the beta parameters, in which each column corresponds to a transition ($1 \rightarrow 2, 1 \rightarrow 3, 1 \rightarrow 4, 2 \rightarrow 1, \dots$), and each row corresponds to a covariate (intercept, distance to center, time since departure). We set to \verb|NA| the columns of unconstrained transition probabilities, and we again fix the intercept of the other columns to a large negative number (here $-100$) to set the corresponding transition probabilities to be virtually zero (i.e.\ impossible transition).

<<fixpar-ses, echo=TRUE, eval=FALSE>>=
fixbeta <- matrix(c(NA,-100,-100,-100,NA,NA,-100,NA,-100,-100,-100,-100,
                    NA,   0,   0,   0, 0, 0,   0, 0,   0,   0,   0,   0,
                     0,   0,   0,   0, 0,NA,   0, 0,   0,   0,   0,   0),
                  nrow=3,byrow=TRUE)
@ 

Biased random walks are used to model the movement in states 1 and 4, with repulsion away from the colony in the outbound trip, and attraction towards the colony in the inbound trip. For that purpose, we include `colony.angle' as a covariate on the angle mean of the von Mises distributions for turning angles in states 1 and 4. 

<<crw-ses-1, echo=TRUE, eval=FALSE>>=
angleFormula <- ~ state1(colony.angle) + state4(colony.angle)
@ 

To specify the direction of the bias (away from or towards the colony), we fix the parameters linking the mean turning angle to the direction of the colony. In state 1, we fix the coefficient to a large negative number (e.g.\ -100) so that the mean direction is away from the colony. In state 4, we fix it to a large positive number (e.g.\ 100) so that the mean direction is towards the colony. Fixing these parameters to large positive (or negative) values effectively removes the correlated random walk component of Eq. \ref{eq:circ} for states 1 and 4. The four other parameters correspond to the angle concentrations and should be estimated (NAs in \verb|fixPar|).
<<crw-ses-2, echo=TRUE, eval=FALSE>>=
fixPar <- list(angle=c(-100,100,NA,NA,NA,NA),beta=fixbeta)
@ 
Because no covariates are specified for the mean angle of state 2 (searching) and state 3 (foraging), these states are reduced to correlated random walks with a mean turning angle of zero (i.e. $\text{atan2}(0,1)=0$ as in Eq. \ref{eq:circ}).

We can now fit the second model with starting parameter values extracted from the simpler model using \verb|getPar0|. In \verb|fitHMM|, we use the arguments \verb|estAngleMean| and \verb|circularAngleMean| to indicate that the angle mean is to be estimated using circular-circular regression.

<<fit-ses, echo=TRUE, eval=FALSE>>=
Par0 <- getPar0(model=m1, nbStates=4, 
                DM=list(angle=list(mean=angleFormula, concentration=~1)), 
                estAngleMean=list(angle=TRUE), 
                circularAngleMean=list(angle=TRUE), formula=formula)

m2 <- fitHMM(data=data, nbStates=4, dist=list(step="gamma",angle="vm"), 
             Par0=list(step=Par0$Par$step, angle=Par0$Par$angle),
             beta0=Par0$beta, fixPar=fixPar, formula=formula, 
             DM=list(angle=list(mean=angleFormula, concentration=~1)), 
             estAngleMean=list(angle=TRUE), circularAngleMean=list(angle=TRUE), 
             stateNames = stateNames)
@ 

Instead of relying entirely on \verb|fixPar| for parameter constraints, an equivalent model for the transition probabilities could be specified using the special function \verb|betaCol| in \verb|formula|:
<<crw-ses-3, echo=TRUE, eval=FALSE>>=
formula <- ~ betaCol1(colony.dist) + betaCol6(time)

fixbeta <- matrix(c(NA,-100,-100,-100,NA,NA,-100,NA,-100,-100,-100,-100,
                    rep(NA,12),
                    rep(NA,12)),
                  nrow=3,byrow=TRUE)

fixPar <- list(angle=c(-100,100,NA,NA,NA,NA),beta=fixbeta)
@
Here \verb|betaCol1(colony.dist)| specifies an effect of distance to colony only on the transition from state 1 to state 2 (which corresponds to the first column of the beta matrix) and \verb|betaCol6(time)| specifies an effect of time since departure only on the transition from state 2 to state 4 (which corresponds to the sixth column of the beta matrix).  When the special function \verb|betaCol| is used, then \verb|fitHMM| automatically fixes the appropriate elements in the second (`colony.dist') and third (`time') rows of the beta matrix to zero (without the user needing to do so manually using \verb|fixPar|). However, note that the first row (corresponding to the intercept terms) must still be manually fixed to achieve the desired constraints on the transition probability matrix.

\subsubsection{Model 3}
In addition to the covariates included in model 2, we add the effect of distance to colony on the step length and turning angle concentration parameters for states 1 and 4. We specify the following formulas:
<<mod3-ses, echo=TRUE, eval=FALSE>>=
distFormula <- ~ state1(colony.dist) + state4(colony.dist)
stepDM <- list(mean=distFormula, sd=distFormula)
angleDM <- list(mean=angleFormula, concentration=distFormula)
@ 

The initial parameters are extracted from model 2, again using the function \verb|getPar0|. Instead of fixing the mean direction of movement like in model 2, we estimate it here as a trade-off between short-term directional persistence and bias toward (or away) from the colony (i.e. a biased correlated random walk as in Eq. \ref{eq:circ}).
<<fit-ses-2, echo=TRUE, eval=FALSE>>=
# remove fixed angle parameters
fixPar <- list(beta=fixbeta)

# get starting parameters from m2
Par0 <- getPar0(model=m2, nbStates=4, 
                DM = list(step=stepDM, angle=angleDM), 
                estAngleMean=list(angle=TRUE), 
                circularAngleMean=list(angle=TRUE), 
                formula=formula)

# the bias is estimated rather than fixed
Par0$Par$angle[c("mean_1:(colony.angle)","mean_4:(colony.angle)")] <- 0

m3 <- fitHMM(data=data, nbStates=4, dist=list(step="gamma",angle="vm"), 
             Par0=list(step=Par0$Par$step, angle=Par0$Par$angle),
             beta0=Par0$beta, fixPar=fixPar, formula=formula, 
             DM = list(step=stepDM, angle=angleDM), 
             estAngleMean=list(angle=TRUE), 
             circularAngleMean=list(angle=TRUE), 
             stateNames = stateNames)
@ 

The three fitted model can be compared with \verb|AIC(m1,m2,m3)|, which overwhelmingly supports model 3. The most likely state sequence is obtained with \verb|viterbi(m3)|. Figure \ref{fig:sesTracks} shows a map of the state-decoded track. The estimated circular-circular regression coefficients for the angle means of state 1 (outbound) and state 4 (inbound) were $\Sexpr{round(sesCIbeta$angle$est[,"mean_1:(colony.angle)"],2)}$ (95\% CI: $\Sexpr{round(sesCIbeta$angle$lower[,"mean_1:(colony.angle)"],2)}-\Sexpr{round(sesCIbeta$angle$upper[,"mean_1:(colony.angle)"],2)}$) and $\Sexpr{round(sesCIbeta$angle$est[,"mean_4:(colony.angle)"],2)}$ (95\% CI: $\Sexpr{round(sesCIbeta$angle$lower[,"mean_4:(colony.angle)"],2)}-\Sexpr{round(sesCIbeta$angle$upper[,"mean_4:(colony.angle)"],2)}$), respectively, thus indicating biased correlated random walks with repulsion away from the colony during outbound movements and attraction towards the colony during inbound movements. The estimated regression coefficients for the step length mean and turning angle concentration parameters for states 1 and 4 suggest that step lengths decreased and turning angles became more concentrated at the mean angle as distance to colony increased (Figure \ref{fig:sesResults})

\begin{figure}[htbp]
  \centering
  \includegraphics[width=\textwidth]{plot_sesResults2.png}
  \caption{The 15 elephant seal tracks, colored by the most likely state sequence.}
  \label{fig:sesTracks}
\end{figure}

\begin{figure}[htbp]
  \centering
  \includegraphics[width=0.49\textwidth]{plot_sesResults001.pdf}
  \includegraphics[width=0.49\textwidth]{plot_sesResults003.pdf}\\
  \includegraphics[width=0.49\textwidth]{plot_sesResults007.pdf}
  \includegraphics[width=0.49\textwidth]{plot_sesResults009.pdf}\\
  \caption{Selected results from the elephant seal example. Panels include estimates and 95\% confidence intervals for the ``outbound'' mean step length parameter (top-left panel), ``inbound'' mean step length parameter (top-right panel), ``outbound'' turning angle concentration parameter (bottom-left panel), and ``inbound'' turning angle concentration parameter (bottom-right panel) as a function of distance to colony (`colony.dist'). Distance to colony has been standardized based on a mean of 2539 km (SD = 1021.3).}
  \label{fig:sesResults}
\end{figure}

\subsection{Group dynamic animal movement}
\label{sec:groupModel}
Here we demonstrate how \verb|momentuHMM| can be used to simulate and fit the group dynamic animal movement model of \cite{LangrockEtAl2014}. In group dynamic models, groups (e.g., herds, packs, schools) are allowed to influence the movement of social individuals. One way to accomplish this is to model individual movements as being attracted to a group ``centroid''. Depending on the system, the centroid could simply be the location of the group center (e.g., the mathematical centroid of the group) or group leader (e.g., alpha wolf) at times $t=1,\dots,T$.  In this sense, the centroid can be considered a dynamic activity center that changes position over time, and these models are not necessarily limited to groups.  For example, the centroid could instead refer to predators, competitors, or human activity (in which case the centroid might be repulsive rather than attractive!). 

Dynamic activity centers can be simulated in \verb|simData| using the \verb|centroids| argument.  Following simulation scenario A of \cite{LangrockEtAl2014}, we first simulate a group centroid as a single-state (i.e. $N=1$) biased correlated random walk relative to the origin:
<<groupModel-centroid, echo=TRUE, eval=FALSE>>=
dist <- list(step="gamma", angle="vm")
nbObs <- 250

Parc <- list(step = c(15,10), 
             angle = c(0.15,log(1)))

DMc <- list(angle=list(mean = ~center1.angle,
                       concentration=~1))

centroidData <- simData(nbStates=1, dist=dist, Par=Parc, DM=DMc,
                        circularAngleMean = list(angle = TRUE),
                        centers = matrix(0,1,2),
                        obsPerAnimal = nbObs)
@
\begin{figure}[htbp]
  \includegraphics[width=0.49\textwidth]{plot_groupExampleCentroid001.pdf}
  \includegraphics[width=0.49\textwidth]{plot_groupExample001.pdf}\\
  \includegraphics[width=0.49\textwidth]{plot_groupExampleResults007.pdf}
  \includegraphics[width=0.49\textwidth]{plot_groupExampleResults008.pdf}\\
  \caption{Selected results from the group dynamic animal movement example. Panels include the simulated centroid path (top-left panel), the simulated paths of 20 individuals where state 1 (``group'') includes biased movements towards the centroid and state 2 (``solitary'') is a correlated random walk independent of the group centroid (top-right panel), and two fitted tracks that are colored by the most likely state sequence (bottom panels).}
  \label{fig:groupResults}
\end{figure}
Now we can use the simulated centroid track (Fig. \ref{fig:groupResults}) as a dynamic activity center and simulate the movement of a group of 20 individuals as a 2-state mixture of a biased random walk (relative to the centroid) and a correlated random walk (independent of the centroid):
<<groupModel-groupData, echo=TRUE, eval=FALSE>>=
nbAnimals <- 20
nbStates <- 2
stateNames <- c("group","solitary")

Par <- list(step = c(30,50,15,25), 
            angle = c(1.e+7,log(2.5),log(5)))

beta <- matrix(c(-2.944439,-1.734601),1,nbStates)

DM <- list(angle=list(mean = ~state1(centroid.angle),
                      concentration = ~1))

# calculate stationary distribution
gamma <- diag(nbStates)
gamma[!gamma] <- exp(beta)
gamma <- t(gamma)
gamma <- gamma/apply(gamma,1,sum)
delta <- solve(diag(nbStates) - t(gamma) + 1, rep(1, nbStates))

# draw random initial locations for each individual
initialPositions <- vector("list")
for (i in 1:nbAnimals) {
  initialPositions[[i]] <- runif(2, -10, 10)
}

# create centroid data frame
cD <- data.frame(x = centroidData$x, y = centroidData$y)

groupData <- simData(nbAnimals=nbAnimals, nbStates=nbStates, dist=dist,
                     Par = Par, beta = beta, delta = delta, DM = DM,
                     circularAngleMean = list(angle = TRUE),
                     centroids = list(centroid = cD),
                     obsPerAnimal = nbObs,
                     initialPosition = initialPositions,
                     states = TRUE, stateNames = stateNames)
@
\noindent Here state 1 (``group'') has biased movements toward the centroid and state 2 (``solitary'') is simply a correlated random walk independent of the group centroid (Fig. \ref{fig:groupResults}).  Note that despite this being a 2-state HMM, the working scale parameters for turning angle (\verb|Par$angle|) only includes 3 parameters (1 for the angle mean and 2 for the concentration parameters).  This is because under the circular-circular regression model, no working parameter is specified\footnote{More accurately, the working parameter for the reference angle is automatically fixed to $\beta_0=1$.} for the reference turning angle of zero (i.e., the component for short-term directional persistence; see Eq. \ref{eq:circ}) and no angular covariates were specified in the model for state 2 (``solitary''). Thus the first parameter corresponds to the working scale parameter of the \verb|centroid.angle| covariate for state 1 (``group''), while the second and third parameters are the working scale parameters for the concentration parameters for states 1 and 2, respectively.  In this case, the angle mean parameter for state 1 was set at a large value (1.e+7) so as to enforce a biased random walk (instead of a biased correlated random walk). In other words, the large value for the biased random walk component ``swamps'' the correlated random walk component in Eq. \ref{eq:circ}.

Finally, we can fit the group dynamic model using \verb|fitHMM|:
<<groupModel-fit, echo=TRUE, eval=FALSE>>=
Par0 <- list(step = c(30,50,15,25), 
             angle = c(0,log(2.5),log(5)))

groupFit <- fitHMM(groupData, nbStates=nbStates, dist=dist, Par=Par0,
                   DM = DM, stationary = TRUE,
                   estAngleMean = list(angle = TRUE),
                   circularAngleMean = list(angle = TRUE),
                   stateNames = stateNames)
@

\subsection{Harbour seals}
\label{sec:harbourSeal}
Here we demonstrate how more complicated parameter constraints can be implemented using the \verb|userBounds| and \verb|workBounds| arguments in \verb|fitHMM| and \verb|MIfitHMM|. This example is based on the harbour seal analysis of \cite{McClintockEtAl2013c}. Using individual-level random effects on probability distribution parameters, \cite{McClintockEtAl2013c} performed a Bayesian analysis of population-level activity budgets for 3-states (``resting'', ``foraging'', and ``transit''). While \verb|momentuHMM| cannot be used to replicate this analysis exactly, we can perform a similar analysis in the absence of individual-level random effects. Here we will focus on several specific parameter constraints, but the full example code can be found in the ``vignettes'' source directory. 

The harbour seal data consist of 17 individuals (10 male, 7 female) and, as in the northern fur example in section \ref{sec:nfs}, both location and dive activity data. The location data were obtained at temporally-irregular intervals, while the dive activity data were obtained at regular 2-hour time steps. We therefore first used \verb|crawlWrap| to fit and predict locations for all 17 tracks at 2-hour time steps and then used \verb|crawlMerge| to merge the predicted locations with the dive activity data. We then fitted several different models assuming: 1) no individual- or sex-level effects on all parameters (i.e., the ``null'' model); 2) sex-level effects on all parameters; and 3) individual-level effects on all parameters.  Based on \verb|fitHMM| fits for the best predicted tracks, AIC overwhelmingly supported the model including individual-level effects on all parameters, but for simplicity we will use the model including no individual- or sex-level effects to demonstrate how the constraints of \cite{McClintockEtAl2013c} can be implemented in \verb|momentuHMM|.  In addition to the lack of individual-level random effects in our example, we also depart from \cite{McClintockEtAl2013c} in our use of zero-inflation parameters to account for steps of length zero (i.e., $l_t=0$) and time steps with no dive activity (i.e., $\omega_t=0$). Unlike \cite{McClintockEtAl2013c}, note that our model 3 also includes individual-level fixed effects on the state transition probabilities.

\cite{McClintockEtAl2013c} fit their 3-state model using more complicated constraints on the probability distribution parameters than any of our previous vignette examples. In addition to relational constraints among the ``resting'', ``foraging'', and ``transit'' states similar to those used in the northern fur seal example (section \ref{sec:nfs}), these constraints included upper bounds for the shape and scale parameters of the Weibull distribution for step length, a minimum value for the ``transit'' concentration parameter of the wrapped Cauchy distribution for turning angle, and bounds on the shape parameters of the beta distribution for dive activity specifically chosen to prevent any ``bathtub'' shaped distributions for the proportion of each time step spent diving below 1.5m. 

Before we demonstrate how to implement these constraints, we should first provide more detail on exactly how the \verb|DM|, \verb|userBounds|, and \verb|workBounds| arguments work together in \verb|fitHMM| (and \verb|MIfitHMM|). While the \verb|DM| argument should now be familiar, we have thus far spent little time discussing the latter two arguments. The \verb|userBounds| argument specifies the lower and upper bounds for the natural scale parameters as a 2-column matrix.  By default all working scale parameters $({\boldsymbol \beta}_\theta)$ are bounded on the real line, but the \verb|workBounds| argument can be used to specify the lower and upper bounds for the working scale parameters as a 2-column matrix.  Specifically, \verb|momentuHMM| calculates natural scale parameters with finite bounds as
\begin{eqnarray}
\label{eq:link}
{\boldsymbol \theta} &=&  \left({\mathbf U}_\theta - {\mathbf L}_\theta \right) g^{-1} \left( {\mathbf X}_\theta  {\boldsymbol \beta}_\theta^* \right) + {\mathbf L}_\theta, \nonumber
%A^{\circ\frac12}
\end{eqnarray}
where %${\mathbf W}_\theta$ is a $K$-vector of \verb|workcons| values, ${\mathbf C}_\theta$ is a $K$-vector of \verb|cons| values specifying the element-wise power to which to raise each of the working scale parameters $({\boldsymbol \beta}_\theta)$, 
${\mathbf L}_\theta$ is the lower bound on the natural scale and ${\mathbf U}_\theta$ is the upper bound on the natural scale.  Note that ${\boldsymbol \beta}_\theta^*={\boldsymbol \beta}_\theta$ under the default values for \verb|workBounds|. For natural scale parameters with finite lower bounds and infinite upper bounds, we have
\begin{eqnarray*}
{\boldsymbol \theta} =  g^{-1} \left( {\mathbf X}_\theta  {\boldsymbol \beta}_\theta^* \right) + {\mathbf L}_\theta.
%A^{\circ\frac12}
\end{eqnarray*}

When \verb|workBounds| is specified, then additional link functions are used on the working scale parameters.  For example, for working scale parameters with finite bounds, we have
\begin{eqnarray*}
{\boldsymbol \beta}_\theta^* &=&  \left({\mathbf U}_{\beta_\theta} - {\mathbf L}_{\beta_\theta} \right) \text{logit}^{-1} \left( {\boldsymbol \beta}_\theta \right) + {\mathbf L}_{\beta_\theta},
\end{eqnarray*}
where ${\mathbf L}_{\beta_\theta}$ is the lower bound on the working scale and ${\mathbf U}_{\beta_\theta}$ is the upper bound on the working scale. When constraining working parameters with finite lower bounds (e.g. zero) and infinite upper bounds, we have
\begin{eqnarray*}
{\boldsymbol \beta}_\theta^* =  \exp \left( {\boldsymbol \beta}_\theta \right) + {\mathbf L}_{\beta_\theta}.
\end{eqnarray*}
When constraining working parameters with infinite lower bounds and finite upper bounds (e.g. zero), we have
\begin{eqnarray*}
{\boldsymbol \beta}_\theta^* =  - \left( \exp \left( -{\boldsymbol \beta}_\theta \right) - {\mathbf U}_{\beta_\theta} \right).
\end{eqnarray*}
Although optimization within \verb|fitHMM| and \verb|MIfitHMM| is always performed on ${\boldsymbol \beta}_\theta$, note that ${\boldsymbol \beta}_\theta^*$ (and a delta method approximation for the variance of this transformation) is returned by \verb|CIbeta|, \verb|MIpool|, and \verb|print| function calls.  

For the Weibull distribution parameters for step length, \cite{McClintockEtAl2013c} constrained the shape parameter $a_s<5$ (to prevent too ``peaked'' distributions) and the scale parameter less than the maximum distance a harbour seal could travel in 2 hours at 2 m/s (i.e., $b_s<14400$m) for $s \in \{1=\text{``resting''},2=\text{``foraging''},3=\text{``transit''} \}$. They also constrained $b_3\ge b_2 \ge b_1$.  This is easily accomplished in \verb|momentuHMM| using the \verb|DM|, \verb|userBounds|, and \verb|workBounds| arguments:
<<spec-hsStep, echo=TRUE, eval=TRUE>>=
nbStates <- 3
stateNames <- c("resting", "foraging", "transit")
dist <- list(step = "weibull", angle = "wrpcauchy", dive = "beta")
stepDM<-matrix(c(1,0,0,0,0,0,0,0,0,
                 0,1,0,0,0,0,0,0,0,
                 0,0,1,0,0,0,0,0,0,
                 0,0,0,1,0,0,0,0,0,
                 0,0,0,1,1,0,0,0,0,
                 0,0,0,1,1,1,0,0,0,
                 0,0,0,0,0,0,1,0,0,
                 0,0,0,0,0,0,0,1,0,
                 0,0,0,0,0,0,0,0,1),nrow=3*nbStates,byrow=TRUE,
        dimnames=list(c(paste0("shape_",1:nbStates),
                        paste0("scale_",1:nbStates),
                        paste0("zeromass_",1:nbStates)),
                      c(paste0("shape_",1:nbStates,":(Intercept)"),
                        "scale:(Intercept)","scale_2","scale_3",
                        paste0("zeromass_",1:nbStates,":(Intercept)"))))
stepworkBounds<-matrix(c(rep(-Inf,4),0,0,rep(-Inf,3),
                   rep(Inf,ncol(stepDM))),ncol(stepDM),2,
                   dimnames=list(colnames(stepDM),c("lower","upper")))
stepBounds<-matrix(c(0,5,
                     0,5,
                     0,5,
                     0,14400,
                     0,14400,
                     0,14400,
                     0,1,
                     0,1,
                     0,1),nrow=3*nbStates,byrow=TRUE,
                    dimnames=list(rownames(stepDM),c("lower","upper")))
@
When included in \verb|workBounds| and \verb|userBounds| for the step length data stream, `stepworkBounds' and `stepBounds' above constrain the parameters $\beta^*_{l,5}>0$ and $\beta^*_{l,6}>0$ such that $14400 \ge b_3 \ge b_2 \ge b_1 \ge 0$:
\begin{eqnarray*}
b_1 &=& \left(14400-0\right) \text{logit}^{-1} \left( \beta_{l,4} \right) + 0 \\
b_2 &=& \left(14400-0\right) \text{logit}^{-1} \left( \beta_{l,4} + \exp(\beta_{l,5}) + 0 \right) + 0 \\
b_3 &=& \left(14400-0\right) \text{logit}^{-1} \left( \beta_{l,4} + \exp(\beta_{l,5}) + 0 + \exp(\beta_{l,6}) + 0\right) + 0.
\end{eqnarray*}

In order to force the ``transit'' state to have strong directional persistence, \cite{McClintockEtAl2013c} constrained the concentration parameter $\rho_3>0.75$. In the absence of relational constraints, \verb|userBounds| could be used to constrain $\rho_3>0.75$.  However, because we also wish to constrain  $\rho_2 \le \rho_3$, we must make use of the \verb|workBounds| argument.  We can constrain $\rho_3 \ge 0.75$, $\rho_2 \le \rho_3$, and $\rho_s \le 0.95$ $(s \in \{ 1,2,3 \})$ using the following combination of \verb|DM|, \verb|userBounds|, and \verb|workBounds| arguments:
<<spec-hsAngle, echo=TRUE, eval=TRUE>>=
angleDM <- matrix(c(1,0,0,
                    0,1,1,
                    0,1,0),nrow=nbStates,byrow=TRUE,
                  dimnames=list(paste0("concentration_",1:nbStates),
                                c("concentration_1:(Intercept)",
                                  "concentration_23:(Intercept)",
                                  "concentration_2")))
angleBounds <- matrix(c(0,0.95,
                        0,0.95,
                        0,0.95),nrow=nbStates,byrow=TRUE,
                      dimnames=list(rownames(angleDM),c("lower","upper")))
transitcons <- boot::logit((0.75 - angleBounds[3,1])
                           /(angleBounds[3,2] - angleBounds[3,1]))
angleworkBounds <- matrix(c(-Inf,transitcons,-Inf,
                            rep(Inf,2),0),ncol(angleDM),2,
                          dimnames=list(colnames(angleDM),c("lower","upper")))
@
\noindent When `angleworkBounds' and `angleBounds' are respectively included in \verb|workBounds| and \verb|userBounds| for the turning angle data stream, this yields
\begin{eqnarray*}
\rho_1 &=& \left(0.95-0\right) \text{logit}^{-1} \left( \beta_{\phi,1} \right) + 0 \\
\rho_2 &=& \left(0.95-0\right) \text{logit}^{-1} \left( \exp(\beta_{\phi,2}) + \Sexpr{round(transitcons,2)} - \left( \exp(-\beta_{\phi,3} ) - 0 \right) \right) + 0 \\
\rho_3 &=& \left(0.95-0\right) \text{logit}^{-1} \left( \exp(\beta_{\phi,2}) + \Sexpr{round(transitcons,2)} \right) + 0.
\end{eqnarray*}
Here we constrained $\rho_s \le 0.95$ to avoid numerical convergence issues that can arise with sparse data sets as $\rho_s \rightarrow 1$. Also note that when using \verb|workBounds| to enforce a specific constraint on the natural scale, \verb|userBounds| should not also include the corresponding natural parameter constraint(s). For example, because we are constraining $\rho_3 \ge 0.75$ with the \verb|workBounds| argument, we did not also include a $0.75$ lower bound for $\rho_3$ in `angleBounds' above.

For the beta distribution of the dive activity data, \cite{McClintockEtAl2013c} constrained the shape1 $(\upsilon_s)$ and shape2 $(\delta_s)$ parameters as follows:
\begin{eqnarray*}
1 \le \upsilon_1 \le \delta_1 \le 10 \\
1 \le \delta_2 \le \upsilon_2 \le 10
\end{eqnarray*}
where $\upsilon_2=\upsilon_3$ and $\delta_2=\delta_3$. These constraints can be imposed using the following combination of \verb|DM|, \verb|userBounds|, and \verb|workBounds| arguments:
<<spec-hsDive, echo=TRUE, eval=FALSE>>=
omegaDM <- matrix(c(1,0,0,0,0,0,
                    0,0,1,1,0,0,
                    0,0,1,1,0,0,
                    1,1,0,0,0,0,
                    0,0,1,0,0,0,
                    0,0,1,0,0,0,
                    0,0,0,0,1,0,
                    0,0,0,0,0,1,
                    0,0,0,0,0,1),nrow=nbStates*3,byrow=TRUE,
                dimnames=list(c(paste0("shape1_",1:nbStates),
                                paste0("shape2_",1:nbStates),
                                paste0("zeromass_",1:nbStates)),
                              c("shape_1:(Intercept)","shape2_1",
                                "shape_2:(Intercept)","shape1_2",
                                "zeromass_1:(Intercept)",
                                "zeromass_23:(Intercept)")))
omegaworkBounds <- matrix(c(-Inf,0,-Inf,0,-Inf,-Inf,
                            rep(Inf,ncol(omegaDM))),ncol(omegaDM),2,
                          dimnames=list(colnames(omegaDM),c("lower","upper")))
omegaBounds <- matrix(c(1,10,
                        1,10,
                        1,10,
                        1,10,
                        1,10,
                        1,10,
                        0,1,
                        0,1,
                        0,1),nrow=nbStates*3,byrow=TRUE,
                      dimnames=list(rownames(omegaDM),c("lower","upper")))
@
Lastly, we wish to impose some biologically-meaningful constraints on the zero-inflation parameters using the \verb|fixPar| argument. Because we would never expect a harbour seal in the ``transit'' state to exhibit a step length of zero, it makes sense to constrain the zero-mass step length parameter for the ``transit'' state to (effectively) zero.  Similarly, we would not expect a harbour seal in the ``foraging'' or ``transit'' states to exhibit no dive activity, and it therefore also makes sense to constrain the zero-mass dive activity parameters for these states to zero.  We can accomplish this using the following \verb|fixPar| argument:
<<spec-hsFixPar, echo=TRUE, eval=FALSE>>=
fixPar <- list(step=c(rep(NA,nbStates*2),NA,NA,boot::logit(1.e-100)),
             omega=c(rep(NA,4),NA,boot::logit(1.e-100)))
@
Putting it all together, we can fit our constrained model assuming no individual- or sex-level effects using \verb|MIfitHMM|:
<<fit-hs, echo=TRUE, eval=FALSE>>=
DM <- list(step = stepDM, angle = angleDM, omega = omegaDM)
userBounds <- list(step = stepBounds, 
                   angle = angleBounds, 
                   omega = omegaBounds)
workBounds <- list(step = stepworkBounds, 
                   angle = angleworkBounds, 
                   omega = omegaworkBounds)
hsFits <- MIfitHMM(crwOut, nSims = 30,
                   nbStates = nbStates, dist = dist, Par0 = Par0,
                   DM = DM, workBounds = workBounds, 
                   userBounds = userBounds, workBounds = workBounds,
                   fixPar = fixPar, stateNames = stateNames)
@

As was mentioned earlier, we found overwhelming AIC support for the individual-level effects model relative to the sex-level effects models and the null model above. While our best supported \verb|momentuHMM| model included individual-level fixed effects, the estimated tracks (Fig. \ref{fig:hsTracks}) and inferences about population-level activity budgets were similar to the individual-level random effects model of \cite{McClintockEtAl2013c}. Estimated activity budgets for the males were $\Sexpr{round(hsActivityBudgets$est[which(hsActivityBudgets$est$sex=="M"),"resting"],2)}$ (95\% CI: $\Sexpr{round(hsActivityBudgets$lower[which(hsActivityBudgets$est$sex=="M"),"resting"],2)}-\Sexpr{round(hsActivityBudgets$upper[which(hsActivityBudgets$est$sex=="M"),"resting"],2)}$) for ``resting'', $\Sexpr{round(hsActivityBudgets$est[which(hsActivityBudgets$est$sex=="M"),"foraging"],2)}$ (95\% CI: $\Sexpr{round(hsActivityBudgets$lower[which(hsActivityBudgets$est$sex=="M"),"foraging"],2)}-\Sexpr{round(hsActivityBudgets$upper[which(hsActivityBudgets$est$sex=="M"),"foraging"],2)}$) for ``foraging'', and $\Sexpr{round(hsActivityBudgets$est[which(hsActivityBudgets$est$sex=="M"),"transit"],2)}$ (95\% CI: $\Sexpr{round(hsActivityBudgets$lower[which(hsActivityBudgets$est$sex=="M"),"transit"],2)}-\Sexpr{round(hsActivityBudgets$upper[which(hsActivityBudgets$est$sex=="M"),"transit"],2)}$) for ``transit''.  Activity budgets for females were $\Sexpr{round(hsActivityBudgets$est[which(hsActivityBudgets$est$sex=="F"),"resting"],2)}$ (95\% CI: $\Sexpr{round(hsActivityBudgets$lower[which(hsActivityBudgets$est$sex=="F"),"resting"],2)}-\Sexpr{round(hsActivityBudgets$upper[which(hsActivityBudgets$est$sex=="F"),"resting"],2)}$) for ``resting'', $\Sexpr{round(hsActivityBudgets$est[which(hsActivityBudgets$est$sex=="F"),"foraging"],2)}$ (95\% CI: $\Sexpr{round(hsActivityBudgets$lower[which(hsActivityBudgets$est$sex=="F"),"foraging"],2)}-\Sexpr{round(hsActivityBudgets$upper[which(hsActivityBudgets$est$sex=="F"),"foraging"],2)}$) for ``foraging'', and $\Sexpr{round(hsActivityBudgets$est[which(hsActivityBudgets$est$sex=="F"),"transit"],2)}$ (95\% CI: $\Sexpr{round(hsActivityBudgets$lower[which(hsActivityBudgets$est$sex=="F"),"transit"],2)}-\Sexpr{round(hsActivityBudgets$upper[which(hsActivityBudgets$est$sex=="F"),"transit"],2)}$) for ``transit''. We found considerable individual variation in the state transition probabilities (Fig. \ref{fig:hsTPM}), and when comparing the estimated activity budgets of our analysis with those of \cite{McClintockEtAl2013c}, we suspect the more noticeable differences between the time spent in the ``resting'' and ``foraging'' states for males is attributable to our having included individual-level effects on the state transition probabilities.  
\begin{figure}[htbp]
  \centering
  \includegraphics[width=\textwidth]{plot_harbourSeal.png}
  \caption{Two harbour seal tracks, colored by the most likely state sequence.}
  \label{fig:hsTracks}
\end{figure}
\begin{figure}[htbp]
  \centering
  \includegraphics[width=\textwidth]{plot_harbourSealResults020}
  \caption{Estimated individual-level state transition probabilities for the harbour seal example.}
  \label{fig:hsTPM}
\end{figure}

\subsection{Northern fulmars}
\label{sec:northernFulmar}
Using Bayesian analysis methods, \cite{PirottaEtAl2018} fit a 6-state biased random walk model to northern fulmar ({\it Fulmarus glacialis}) tracks in northern Scotland, UK.  These states included biased movements relative to a colony in Eynhallow ($59.12^\circ$ N, $3.1^\circ$ W) and fishing vessels that frequently work in the area.  \cite{PirottaEtAl2018} framed their model as having two latent state process.  Under the first state process, the direction of movement could be biased away from the colony (``sea''), towards the nearest fishing vessel (``boat''), or towards the colony (``colony''). Under the second state process, the movement mode could either be fast and directionally-persistent (``transit'') or area-restricted search (``ARS'').  Thus the six states are ``sea ARS'', ``sea transit'', ``boat ARS'', ``boat transit'', ``colony ARS'' and ``colony transit''. \cite{PirottaEtAl2018} also allowed the distance to the nearest fishing vessel and time since leaving the colony to affect state transitions to the ``boat'' and ``colony'' states.  Here we demonstrate how a very similar (but not identical) model can be implemented in \verb|momentuHMM|.

The data are provided in a Dryad repository, but these will require some additional formatting and preparation.  We first load the raw data, create the required ``ID'' column based on individual trips, convert time stamps to class \verb|POSIXct|, and project the northern fulmar (``Longitude'', ``Latitude'') and nearest fishing vessel (``Boat\_Longitude'', ``Boat\_Latitude'') locations using the \verb|sp| package:
<<prep-fulmar-1, echo=TRUE, eval=FALSE>>=
library(sp)

# load data provided by Pirotta et al
fulmarURL <- "https://datadryad.org/bitstream/handle/10255/dryad.174482/"
fileName <- "Fulmar_trackingData.csv?sequence=1"
raw_data <- read.csv(url(paste0(fulmarURL,fileName)),
                     stringsAsFactors = FALSE)

raw_data$ID <- raw_data$tripID
raw_data$Date <- as.POSIXct(raw_data$Date,tz="UTC",
                            format="%d/%m/%Y %H:%M")

# project data
oldProj <- CRS("+proj=longlat +datum=WGS84")
newProj <- CRS("+init=epsg:27700")
coordinates(raw_data) <- c("Longitude","Latitude")
proj4string(raw_data) <- oldProj
raw_data <- as.data.frame(spTransform(raw_data, newProj))

coordinates(raw_data) <- c("Boat_Longitude","Boat_Latitude")
proj4string(raw_data) <- oldProj
raw_data <- as.data.frame(spTransform(raw_data, newProj))
@
For movements away from the colony (``sea ARS'' and ``sea transit''), \cite{PirottaEtAl2018} included bias in the direction of the farthest location from the colony for a given trip. We can use the \verb|centers| argument of \verb|prepData| to identify these locations  (``max\_dist'') and then calculate the expected angle for the ``sea'' states (``sea.angle'') using \verb|momentuHMM:::distAngle|:
<<prep-fulmar-2, echo=TRUE, eval=FALSE>>=
# use prepData to calculate colony distance covariate ('sea.angle')
colony <- data.frame(x = -3.1, y = 59.12) 
coordinates(colony) <- c("x", "y")
proj4string(colony) <- oldProj
colony <- as.matrix(as.data.frame(spTransform(colony, newProj)))
rownames(colony) <- "colony"
colony_dist <- prepData(raw_data, coordNames = c("Longitude","Latitude"),
                        centers = colony)

# calculate "sea" mean angle covariate
sea.angle <- NULL
for(id in unique(colony_dist$ID)) {
  idat <- subset(colony_dist,ID==id)
  nbSubObs <- length(which(colony_dist$ID==id))
  max_dist <- as.numeric(idat[which.max(idat$colony.dist),c("x","y")])
  max_angle <- momentuHMM:::distAngle(colony,colony,max_dist)[2]
  sea.angle <- c(sea.angle, rep(max_angle,nbSubObs))
}
raw_data$sea.angle <- sea.angle
@

Next we calculate the time since leaving colony covariate (``time''):
<<prep-fulmar-3, echo=TRUE, eval=FALSE>>=
# calculate time since left colony covariate ('time')
time <- aInd <- NULL
for(id in unique(raw_data$ID)) {
  idInd <- which(raw_data$ID==id)
  aInd <- c(aInd,idInd[1])
  nbSubObs <- length(idInd)
  time <- c(time, (1:nbSubObs)/nbSubObs)
}
raw_data$time <- time
@

To complete our data preparation, we convert the nearest fishing vessel data to the \verb|centriods| argument format and use \verb|prepData| to calculate step lengths, turn angles, and our ``sea'', ``boat'', and ``colony'' covariates:
<<prep-fulmar-4, echo=TRUE, eval=FALSE>>=
# get boat data into centroids argument format
boat_data <- list(boat=data.frame(Date = raw_data$Date,
                                  x = raw_data$Boat_Longitude,
                                  y = raw_data$Boat_Latitude))

# format and merge all data and covariates for analysis
fulmar_data <- prepData(raw_data, coordNames = c("Longitude","Latitude"),
                        centers = colony,
                        centroids = boat_data,
                        covNames = "time",
                        angleCovs = "sea.angle")

# momentuHMM doesn't like data streams and covariates to have same name, 
# so create identical data column with different name
fulmar_data$d <- fulmar_data$boat.dist 

# standarize boat.dist covariate
fulmar_data$boat.dist <- scale(fulmar_data$boat.dist)
@
Note that we use the \verb|centers| argument for the colony (because its location is static) and the \verb|centroids| argument for the nearest fishing vessel (because its location is dynamic).  

Now that we've formatted the data, we're ready to specify the 6-state HMM. Using 10 min time steps, \cite{PirottaEtAl2018} included three data streams in their model: step length (``step''), turn angle (``angle''), and distance to nearest boat (``d'').  These were respectively modelled using Weibull, wrapped Cauchy, and log-normal distributions:
<<spec-fulmar-1, echo=TRUE, eval=FALSE>>=
nbStates <- 6

stateNames <- c("seaARS", "seaTr", 
                "boatARS", "boatTr", 
                "colonyARS", "colonyTr")

dist <- list(step = "weibull",
             angle = "wrpcauchy",
             d = "lnorm")

@

Similar to the harbour seal example (section \ref{sec:harbourSeal}), \cite{PirottaEtAl2018} used relational parameter constraints that can be specified in \verb|momentuHMM| using pseudo-design matrices:
<<spec-fulmar-2, echo=TRUE, eval=FALSE>>=
# specify data stream probability distribution parameter constraints
stepDM <- matrix(c(1,0,0,0,
                   0,1,0,0,
                   1,0,0,0,
                   0,1,0,0,
                   1,0,0,0,
                   0,1,0,0,
                   0,0,1,0,
                   0,0,1,1,
                   0,0,1,0,
                   0,0,1,1,
                   0,0,1,0,
                   0,0,1,1),2*nbStates,4,byrow=TRUE,
                 dimnames=list(c(paste0("shape_",1:nbStates),
                                 paste0("scale_",1:nbStates)),
                               c("shape:ARS","shape:Tr",
                                 "scale:(Intercept)","scale:Tr")))

# constrain scale parameters such that Tr > ARS 
stepworkBounds <- matrix(c(-Inf,Inf,
                           -Inf,Inf,
                           -Inf,Inf,
                           0,Inf),ncol(stepDM),2,byrow=TRUE,
                         dimnames=list(colnames(stepDM),c("lower","upper")))

# include trip-level effects on angle mean concentration parameter
nbTrips <- length(unique(fulmar_data$ID))
angleDM <- matrix(c("sea.angle",0,0,0,0,rep(0,2*nbTrips),
                    "sea.angle",0,0,0,0,rep(0,2*nbTrips),
                    0,"boat.angle",0,0,0,rep(0,2*nbTrips),
                    0,"boat.angle",0,0,0,rep(0,2*nbTrips),
                    0,0,"colony.angle",0,0,rep(0,2*nbTrips),
                    0,0,"colony.angle",0,0,rep(0,2*nbTrips),
                    0,0,0,1,0,paste0("ID",1:nbTrips),rep(0,nbTrips),         
                    0,0,0,1,1,paste0("ID",1:nbTrips),paste0("ID",1:nbTrips),
                    0,0,0,1,0,rep(0,2*nbTrips),
                    0,0,0,1,1,rep(0,2*nbTrips),
                    0,0,0,1,0,rep(0,2*nbTrips),
                    0,0,0,1,1,rep(0,2*nbTrips)),2*nbStates,3+2+2*nbTrips,byrow=TRUE,
                  dimnames=list(c(paste0("mean_",1:nbStates),
                                  paste0("concentration_",1:nbStates)),
                                c("mean:sea","mean:boat","mean:colony",
                                  "concentration:(Intercept)","concentration:Tr",
                                  paste0("concentration:ID",1:nbTrips,":(Intercept)"),
                                  paste0("concentration:ID",1:nbTrips,":Tr"))))

# constrain concentration parameters such that Tr > ARS 
angleworkBounds <- matrix(c(-Inf,Inf,
                            -Inf,Inf,
                            -Inf,Inf,
                            -Inf,Inf,
                            0,Inf,
                            rep(c(-Inf,Inf),nbTrips),
                            rep(c(0,Inf),nbTrips)),ncol(angleDM),2,byrow=TRUE,
                         dimnames=list(colnames(angleDM),c("lower","upper")))

dDM <- matrix(c(1,1,0,0,
                1,1,0,0,
                1,0,0,0,
                1,0,0,0,
                1,1,0,0,
                1,1,0,0,
                0,0,1,1,
                0,0,1,1,
                0,0,1,0,
                0,0,1,0,
                0,0,1,1,
                0,0,1,1),2*nbStates,4,byrow=TRUE,
              dimnames=list(c(paste0("location_",1:nbStates),
                              paste0("scale_",1:nbStates)),
                            c("location:(Intercept)","location:noboat",
                              "scale:(Intercept)","scale:noboat")))

# constrain location and scale parameters such that sea and colony > boat
dworkBounds <- matrix(c(-Inf,Inf,
                        0,Inf,
                        -Inf,Inf,
                        0,Inf),ncol(dDM),2,byrow=TRUE,
                      dimnames=list(colnames(dDM),c("lower","upper")))

DM <- list(step = stepDM, angle = angleDM, d = dDM)

workBounds <- list(step = stepworkBounds, 
                   angle = angleworkBounds, 
                   d = dworkBounds)
@
To complete our model specification, we use the \verb|toState| special function to model transitions to the ``boat'' and ``colony'' states as a function of distance to nearest fishing vessel (``boat.dist'') and time since leaving colony (``time''), respectively.  Following \cite{PirottaEtAl2018}, we will use the \verb|knownStates| argument to fix the initial state to ``sea transit''.  We will also use the \verb|fixPar| argument to fix the initial state probabilities (because we are assuming these are known) and, as in section \ref{sec:groupModel}, constrain the model to a biased random walk by fixing the mean angle working scale parameters to a large positive value:
<<spec-fulmar-3, echo=TRUE, eval=FALSE>>=
# state transition formula similar to Pirotta et al
formula <- ~ toState3(boat.dist) + toState4(boat.dist) + 
             toState5(time) + toState6(time)

# specify knownStates
# Priotta et al assumed all animals start in state 2 ('seaTr')
knownStates <- rep(NA,nrow(fulmar_data))
knownStates[aInd] <- 2 

# fix delta_2 = 1 because assuming initial state is known for each track
fixPar <- list(delta=c(100,rep(0,nbStates-2)))
fixPar$delta <- exp(c(0,fixPar$delta))/sum(exp(c(0,fixPar$delta)))
# Constrain model to BRW (instead of BCRW)
fixPar$angle <- c(rep(1.e+7, 3), rep(NA, 2+2*nbTrips)) 
@
\noindent Lastly, we used the \verb|betaCons| argument to impose similar constraints as \cite{PirottaEtAl2018} for the transition probability parameters. We accomplish this by using \verb|betaCons| to set the transition probability working parameter intercept terms equal among the three ``ARS'' states and the three ``transit'' states. We also use \verb|betaCons| to constrain the effects of `boat.dist' and `time' to be identical for each of the two movement modes (``ARS'' and ``transit'') within each of the three biased movement states (``sea'', ``boat'', and ``colony''). \verb|betaCons| must be a matrix of the same dimension as \verb|beta0| and be composed of integers, where each beta working parameter is sequentially indexed in a column-wise fashion. Equality constraints can then be incorporated by having parameters share the same index. In this example we have:
<<spec-fulmar-4, echo=FALSE, eval=TRUE>>=
nbStates <- 6
betaCons <- matrix(1:(3*nbStates*(nbStates-1)),3,nbStates*(nbStates-1),
                   dimnames=list(c("(Intercept)","boat.dist","time"),
                                 paste(rep(1:nbStates,each=nbStates),"->",rep(1:nbStates,nbStates))[(1:(nbStates*nbStates))[-diag(matrix(1:(nbStates*nbStates),nbStates,nbStates))]]))
betaCons["(Intercept)",c("1 -> 4","1 -> 6","3 -> 2","3 -> 4","3 -> 6","5 -> 2","5 -> 4","5 -> 6")] <- betaCons["(Intercept)","1 -> 2"] # constrain ARS -> Tr intercept
betaCons["(Intercept)",c("1 -> 5","3 -> 1","3 -> 5","5 -> 1","5 -> 3")] <- betaCons["(Intercept)","1 -> 3"] # constrain ARS -> ARS intercept
betaCons["(Intercept)",c("2 -> 3","2 -> 5","4 -> 1","4 -> 3","4 -> 5","6 -> 1","6 -> 3","6 -> 5")] <- betaCons["(Intercept)","2 -> 1"] # constrain Tr -> ARS intercept
betaCons["(Intercept)",c("2 -> 6","4 -> 2","4 -> 6","6 -> 2","6 -> 4")] <- betaCons["(Intercept)","2 -> 4"] # constrain Tr -> Tr intercept
betaCons["boat.dist",c("1 -> 4","2 -> 3","2 -> 4")] <- betaCons["boat.dist","1 -> 3"] # constrain boat.dist 1 -> 3 = 1 -> 4 = 2 -> 3 = 2 -> 4
betaCons["boat.dist","4 -> 3"] <- betaCons["boat.dist","3 -> 4"] # constrain boat.dist 3 -> 4 = 4 -> 3
betaCons["boat.dist",c("5 -> 4","6 -> 3","6 -> 4")] <- betaCons["boat.dist","5 -> 3"] # constrain boat.dist 5 -> 3 = 5 -> 4 = 6 -> 3 = 6 -> 4
betaCons["time",c("1 -> 6","2 -> 5","2 -> 6")] <- betaCons["time","1 -> 5"] # constrain time 1 -> 5 = 1 -> 6 = 2 -> 5 = 2 -> 6
betaCons["time",c("3 -> 6", "4 -> 5", "4 -> 6")] <- betaCons["time","3 -> 5"] # constrain time 3 -> 5 = 3 -> 6 = 4 -> 5 = 4 -> 6
betaCons["time","6 -> 5"] <- betaCons["time","5 -> 6"] # constrain time 5 -> 6 = 6 -> 5

@
<<spec-fulmar-5, echo=TRUE, eval=TRUE, size='small'>>=
betaCons
@
\noindent Again, \verb|betaCons| constrains any of the transition probability matrix working parameters with the same index to be equal to one another. For example, all of the intercept terms indexed by a `1' ($1 \rightarrow 2$, $1 \rightarrow 4$, $1 \rightarrow 6$, $3 \rightarrow 2$, $3 \rightarrow 4$, $3 \rightarrow 6$, $5 \rightarrow 2$, $5 \rightarrow 4$, and $5 \rightarrow 6$) are equal, and these terms correspond to the transitions from the ``ARS'' movement mode (states 1, 3, and 5) to the ``transit'' movement mode (states 2, 4, and 6). Similarly, all of the intercept terms indexed by a `16' are equal and correspond to the transitions from the ``transit'' movement mode to the ``ARS'' movement mode. For further details on the \verb|betaCons| argument, see the \verb|fitHMM| help file and the northern fulmar example code in the ``vignettes'' source directory.

Now we are ready to fit our 6-state HMM:
<<fit-fulmar-1, echo=TRUE, eval=FALSE>>=
m2 <- fitHMM(fulmar_data, nbStates, dist,
             Par0 = Par0$Par, beta0 = Par0$beta0,
             formula = formula,
             estAngleMean = list(angle = TRUE), 
             circularAngleMean = list(angle = TRUE),
             DM = DM, workBounds = workBounds, betaCons = betaCons,
             fixPar = fixPar, knownStates = knownStates,
             stateNames = stateNames) 
@
\noindent With decent starting values, this model required about \Sexpr{fulmarElapsedTime} min to fit on a standard desktop computer (macOS El Capitan, 2.8 GHz Intel Core i7, 16 GB RAM). For comparison, \cite{PirottaEtAl2018} required about 18 hr to fit their Bayesian model using MCMC (2.9 GHz Intel Core i7, 16 GB RAM). 

We can compare the estimated activity budgets with those of \cite{PirottaEtAl2018} using the \verb|timeInStates| function:
<<results-fulmar-1a, echo=TRUE, eval=FALSE>>=
timeInStates(m2)
@
<<results-fulmar-1b, echo=FALSE>>=
timeIn1
@
\noindent \cite{PirottaEtAl2018} estimated 0.28 (``seaARS''), 0.20 (``seaTr''), 0.18 (``boatARS''), 0.06 (``boatTr''), 0.12 (``colonyARS''), and 0.16 (``colonyTr''). While these are very similar, there a handful of state assignments that differ between the analyses. These differences could be attributable to several factors, including: 1) the use of informative priors in the Bayesian analysis of \cite{PirottaEtAl2018}; 2) our use of fixed trip-level effects on the ``sea'' state turn angle concentration parameters; and 3) \cite{PirottaEtAl2018} assumed the state transition probability covariates (``boat.dist'' and ``time'') only affected the movement direction states (``sea'', ``boat'', ``colony''), but in our \verb|momentuHMM| implementation the covariates can affect state transitions for both the movement direction (``sea'', ``boat'', ``colony'') and the movement mode (``ARS'', ``transit'').

We can also examine activity budgets by individual bird (which are indexed in the raw data ``birdID'' column), where it is clear that the first 3 individuals tended to spend a larger proportion of their foraging trips in the ``boat'' states:
<<results-fulmar-2a, echo=TRUE, eval=FALSE>>=
timeInStates(m2, by = "birdID")
@
<<results-fulmar-2b, echo=FALSE>>=
timeIn2
@

Finally, we can create a plot similar to \cite{PirottaEtAl2018} using the \verb|plotSat| function (Figure \ref{fig:northernFulmar}):
<<results-fulmar-3, echo=TRUE, eval=FALSE>>=
plotSat(m2, zoom = 7, shape = c(17,1,17,1,17,1), size = 2, 
        col = rep(c("#E69F00", "#56B4E9", "#009E73"), each = 2),
        stateNames = c("sea ARS", "sea Transit",
                       "boat ARS", "boat Transit",
                       "colony ARS", "colony Transit"),
        projargs = newProj, ask = FALSE)
@
\begin{figure}[htbp]
  \centering
  \includegraphics[width=\textwidth]{plot_northernFulmarExample.png}
  \caption{Seven northern fulmar tracks, colored by the most likely state sequence.}
  \label{fig:northernFulmar}
\end{figure}

\section{Discussion}
Here we have introduced version \Sexpr{installed.packages()["momentuHMM","Version"]} of the R package \verb|momentuHMM| and demonstrated some of its capabilities for conducting multivariate HMM analyses with animal location, auxiliary biotelemetry, and environmental data. The package allows for fitting (and simulating from) a suite of biased and correlated random walk movement process models \citep[e.g.][]{McClintockEtAl2012}, can be used for an unlimited number of data streams and latent behavior states, includes multiple imputation methods to account for measurement error, temporal irregularity, and other forms of missing data that would otherwise be prohibitive to maximum likelihood analysis, and integrates seamlessly with rasters to facilitate spatio-temporal covariate modelling. Because the package incorporates biased random walks, it can also be used to implement group dynamic models \cite[e.g.][]{LangrockEtAl2014}. The package therefore greatly expands on available software and facilitates the incorporation of more ecological and behavioral realism for hypothesis-driven analyses of animal movement that account for many of the challenges commonly associated with telemetry data. While many of the features of \verb|momentuHMM| were motivated by animal movement data, we note that the package is not limited to location data and can be used for analyzing any type of data that is amenable to (multivariate) HMMs.

Model fitting in \verb|momentuHMM| is relatively fast because the forward algorithm (Eq. \ref{eq:HMMlike}) is coded in C++. Because multiple imputations are completely parallelizable, with sufficient processing power computation times for analyses that account for measurement error, temporal irregularity, or other forms of missing data need not be longer than that required to fit a single HMM.  However, computation times will necessarily be longer as the number of states and/or parameters increase. For example, \verb|momentuHMM| required about 1 hr to fit a single HMM with $N=6$ states, seven data streams, and $T=7414$ time steps \citep{McClintock2017}. 

As in any maximum likelihood analysis based on numerical optimization, computation times will also depend on the starting values (\verb|Par0| and \verb|beta0|). Specifying ``good'' starting values is arguably the most challenging aspect of model fitting in \verb|momentuHMM|, particularly for the working scale coefficients when using covariates. The \verb|getPar|, \verb|getPar0|, \verb|getParDM|, and \verb|checkPar0| functions are designed to help with the specification of starting values, and the \verb|retryFits| argument in \verb|crawlWrap|, \verb|fitHMM|, and \verb|MIfitHMM| will re-optimize based on random perturbations of the parameters to help explore the likelihood surface and diagnose convergence to local maxima. Optimization for the circular-linear regression link function (\verb|tan(mean/2)|; see Table \ref{tab:pdfs}) in particular can be prone to local minima, so users are encouraged to explore a range of starting values when fitting these models.

While \verb|momentuHMM| includes functions for drawing realizations of the position process based on the CTCRW model of \cite{JohnsonEtAl2008}, this is but one of many methods for performing the first stage of multiple imputation. Realizations of the position process from any movement model that accounts for measurement error and/or temporal irregularity \citep[e.g.][]{CalabreseEtAl2016,GurarieEtAl2017} could be passed to \verb|MIfitHMM| for HMM-type analyses in the second stage. Multiple imputation methods also need not be limited to these telemetry error scenarios. For example, conventional missing data could also be imputed using standard techniques \citep{RubinSchenker1986}, thereby allowing the investigation of non-random mechanisms for missingness that can be problematic if left unaccounted for in HMMs.

There remain many potential avenues for refining and extending the capabilities of \verb|momentuHMM|. Computation times could likely be improved by further optimizing the R and C++ code for speed. Notable extensions include hidden semi-Markov models and random effects on data stream probability distribution and state transition probabilitiy parameters \citep{ZucchiniEtAl2016}. We would also like to incorporate additional parameters for change-point thresholds and the locations of activity centers instead of requiring that they be pre-specified (and potentially compared using AIC or other model selection criteria) as in grey seal example. Lastly, it is relatively straightforward to add additional probability distributions, and we are pleased to do so upon request. Practitioners interested in additional features for \verb|momentuHMM| are encouraged to contact the authors.

\section*{Acknowledgments} 
We are grateful to R. Scott, B. Godley, M. Godfrey, J. Sudre, and North Carolina Aquariums for providing the data used in our turtle example. The findings and conclusions in the manuscript are those of the author(s) and do not necessarily represent the views of the National Marine Fisheries Service, NOAA. Any use of trade, product, or firm names does not imply an endorsement by the US Government.

\bibliographystyle{mee}
\bibliography{master}

\clearpage

\end{document}
