Data preparation

The data set is included in moveHMM by default, which includes the first three tracks from the study of Michelot, Langrock, and Patterson (2016). This data set is already in the format expected by moveHMM, with columns:

• ID: an identifier for the track/individual;

• x: Easting or longitude coordinate;

• y: Northing or latitude coordinate.

These three variables are required, and should have these names. The other two columns, slope and temp, are covariates that will be included in the analysis later.

head(haggis_data)

  ID         x         y     slope      temp
1  1  0.000000  0.000000 25.957002 10.344959
2  1 -1.068761 -0.194650 18.606632  8.352531
3  1 -6.152549  2.051343 16.524004 13.529650
4  1 -6.703983  3.338480  9.154917 10.951095
5  1 -6.541667  3.553843  5.547686 11.243328
6  1 -7.160298  1.960377  8.129402 13.187280

The standard HMMs used in movement ecology, and implemented in this package, model step lengths and turning angles. The first step of the analysis is therefore to derive those variables from the observed locations. This can be done using the function prepData; here, we specify type = "UTM" to indicate that the locations are already projected (rather than longitude-latitude locations).

data <- prepData(haggis_data, type = "UTM")

head(data)

  ID      step      angle         x         y     slope      temp
1  1 1.0863417         NA  0.000000  0.000000 25.957002 10.344959
2  1 5.5578218 -0.5961622 -1.068761 -0.194650 18.606632  8.352531
3  1 1.4002860 -0.7500230 -6.152549  2.051343 16.524004 13.529650
4  1 0.2696813 -1.0506197 -6.703983  3.338480  9.154917 10.951095
5  1 1.7093394 -2.8660552 -6.541667  3.553843  5.547686 11.243328
6  1 1.1529149  2.3676683 -7.160298  1.960377  8.129402 13.187280

The output data frame has two new columns: step (step lengths), and angle (turning angle).

Model fitting

In moveHMM, the models are fitted using the function fitHMM. It implements maximum likelihood estimation, using the numerical optimiser nlm. The likelihood function captures how plausible the data are, given a set of parameters, and nlm explores the space of parameter values to find the one that maximise the likelihood. For this reason, it is necessary to choose initial parameter values, from where nlm can start its exploration, and these are specified using the stepPar0 and anglePar0 arguments in fitHMM. The choice of starting values can affect the performance of the optimiser, and in some cases its ability to find the maximum likelihood estimates. Some guidance on how to choose those values is provided in another vignette, called “A short guide to choosing initial parameter values for the estimation in moveHMM”. In summary, one approach is to look at the empirical distributions of step lengths and turning angles, to make an educated guess about what parameter values are plausible.

In the following, we consider a 2-state HMM where, in each state, a gamma distribution models step length, and a von Mises distribution models turning angle. In moveHMM, the gamma distribution is specified in terms of a mean and standard deviation, and we pick some initial values based on the histogram of observed step lengths. The von Mises distribution has a mean and a concentration parameter, and reasonable starting values can be chosen based on a histogram of observed turning angles. Please refer to the dedicated vignette for more guidance about selecting initial parameters.

hist(data$step, xlab = "step length")
hist(data$angle, breaks = seq(-pi, pi, length = 15), xlab = "turning angle")

stepPar0 <- c(1, 5, 1, 5)
anglePar0 <- c(pi, 0, 0.3, 5)

We can fit the model using fitHMM, with arguments nbStates (number of states), stepPar0 (initial step length parameters), and anglePar0 (initial turning angle parameters). By default, the function uses the gamma and von Mises distributions, so we don’t need to specify them here, but the arguments stepDist and angleDist could be included to use different distributions.

mod1 <- fitHMM(data = data, nbStates = 2, 
               stepPar0 = stepPar0, anglePar0 = anglePar0)

mod1

Value of the maximum log-likelihood: -3462.487 

Step length parameters:
----------------------
       state 1  state 2
mean 0.9977756 5.009977
sd   0.4951516 3.044708

Turning angle parameters:
------------------------
                state 1    state 2
mean          -3.109852 -0.3082705
concentration  1.035834  8.7682404

Regression coeffs for the transition probabilities:
--------------------------------------------------
             1 -> 2    2 -> 1
intercept -1.736619 -2.002033

Transition probability matrix:
-----------------------------
          [,1]      [,2]
[1,] 0.8502571 0.1497429
[2,] 0.1189896 0.8810104

Initial distribution:
--------------------
[1] 0.06370057 0.93629943

The object returned by fitHMM is a list that contains, among other things, estimates of all model parameters. Those are automatically shown in a convenient layout when the model object is printed, but they can also be accessed as mod1$mle.

Note that the parameter estimates are pooled across individuals; i.e., a common model is fitted to all tracks in the data set. This is based on the assumption that all animals follow similar movement patterns, and parameters can then be viewed as an average across individuals. In cases where this assumption doesn’t hold, one option is to include individual-specific covariates on the transition probabilities, to capture some of the inter-individual heterogeneity. (See section on covariates below.)

Model visualisation

To understand and interpret a fitted model, it is often helpful to visualise it. The plot function can be called on the model object to generate a few plots, including:

• histograms of the step lengths, overlaid with estimated distributions in each state. This is helpful to check whether the fitted distributions match the data well, and to interpret each state (e.g., does state 1 capture long or short step lengths? how variable are step lengths in state 2?).

• histograms of the turning angles, overlaid with estimated distributions in each state.

• plots of the movement tracks, coloured by estimated state. This provides useful information about where an when animals followed each movement type, what movement patterns were captured in each state, etc.

plot(mod1, ask = FALSE, animals = 1)

Decoding states sequence... DONE

The first state captures short step lengths (mostly between 0 and 2 km) and undirected movement (flat distribution of turning angles), whereas the second state captures long step lengths (between 0 and 15 km) and directed movement (peak of turning angles at zero). In practice, these two states might therefore be used as proxies for two different behavioural states.

State process inference

As shows in the plots above, it is possible to estimate a state for each observation, which is often of great applied interest. The most common way to do this is to use the Viterbi algorithm. The function viterbi returns a sequence of states (i.e., vector of numbers between 1 and nbStates), and these can for example be used to generate plots of the tracks, or of the step lengths, coloured by the most likely state sequence.

# Add most likely state sequence to data
data$state <- factor(viterbi(mod1))

# Plot tracks coloured by state
ggplot(data, aes(x, y, col = state, group = ID)) + 
    geom_path() +
    coord_equal()
# Plot step lengths coloured by state
ggplot(data, aes(x = 1:nrow(data), y = step, col = state, group = ID)) + 
    geom_point(size = 0.3)

An alternative to the Viterbi algorithm is to compute state probabilities, i.e., the probability of being in each state for each observation in the data. State probabilities provide more detailed information about the uncertainty on the state allocation, and they can be computed with the function stateProbs. The output is a matrix with one column for each state and one row for each observation.

sp <- stateProbs(m = mod1)
head(sp)

             [,1]         [,2]
[1,] 9.588326e-02 9.041167e-01
[2,] 1.215803e-06 9.999988e-01
[3,] 3.145690e-01 6.854310e-01
[4,] 9.586125e-01 4.138745e-02
[5,] 1.000000e+00 1.260703e-08
[6,] 1.000000e+00 9.504981e-10

# Add prob of state 1 to data set
data$sp1 <- sp[,1]
ggplot(data, aes(x, y, col = sp1, group = ID)) +
    geom_path() +
    coord_equal() +
    labs(col = "Pr(S = 1)")

Covariates

HMMs have been extended to allow for covariate effects on the transition probabilities. This approach is often interesting to answer ecological questions of the form “does [some covariate] have an effect on the animal’s behaviour?”.

We model the haggises’ transition probabilities as functions of two covariates: temperature, and slope. We expect a non-linear effect of slope, so we include a quadratic term for that covariate. The model formula can be specified using the usual R syntax, and passed through the argument formula in fitHMM. The other arguments are unchanged.

mod2 <- fitHMM(data = data, nbStates = 2, 
               stepPar0 = stepPar0, anglePar0 = anglePar0, 
               formula = ~ temp + slope + I(slope^2))

Plotting the model now also returns plots of the transition probabilities over a grid of values of each covariate (with 95% pointwise confidence intervals if plotCI = TRUE). In this example, it looks like there is no clear effect of temperature, but a clear effect of slope. The probability of remaining in state 2 (or of switching from state 1 to state 2) is highest when the slope is between 10 and 30 degrees, suggesting that wild haggises move faster in that range of slopes. (See Michelot, Langrock, and Patterson (2016) for more details about the interpretation, but it has to do with their morphology.)

plot(mod2, ask = FALSE, plotTracks = FALSE, plotCI = TRUE)

Decoding states sequence... DONE

Another useful output to interpret covariate effects is the stationary state probabilities. These capture the probability of being in each state over the long term, and they can be derived for a grid of covariate values to assess how animals’ activity budget depends on a covariate. The probabilities can be computed with stationary, and plotted with plotStationary. The stationary state probabilities make it even clearer that haggises tend to spend time in state 2 when the slope is between 10 and 30, but are likely to be in state 1 for small or large slope values.

plotStationary(mod2, plotCI = TRUE)

Models can be compared using the AIC function, for example to decide whether or not to keep a covariate. Here, the model with covariate is strongly preferred by AIC.

AIC(mod1, mod2)

  Model      AIC
1  mod2 6814.784
2  mod1 6946.975

Model checking

A chosen HMM formulation rests on many assumptions, including the number of states, the types of distributions used for step length and turning angle, and the dependence structure of the data-generating process. Visual inspection of the model outputs (e.g., fitted distributions) is a first step to check whether it seems to capture features in the data well. Pseudo-residuals are another tool to assess goodness-of-fit in HMMs. Similarly to residuals in linear models, they should be independent and follow a standard normal distribution if all model assumptions were satisfied. Deviations from these patterns suggest lack of fit.

Pseudo-residuals can be computed with pseudoRes, which returns a vector for step lengths and one for turning angles. The function plotPR creates three plots of the pseudo-residuals for each data variable: a time series plot, a quantile-quantile (QQ) plot against the standard normal distribution, and an autocorrelation function (ACF) plot.

plotPR(mod2)

Computing pseudo-residuals... DONE

In the QQ plot, deviations from the diagonal line suggest that the estimated distributions do not capture the empirical distribution of the corresponding variable. Here, the points are aligned on the diagonal, indicating that the fit is good. The ACF plots can be used to detect residual autocorrelation (bars that stretch beyond the confidence band around zero). In this case, the ACF is virtually zero, suggesting that the dependence assumptions of the model are satisfied.

Custom plots

The function getPlotData creates data frames in a format that is convenient for creating custom plots of the model, including state-dependent distributions (type = "dist"), transition probabilities (type = "tpm"), and stationary state probabilities (type = "stat"). The option format specifies whether the data frame should be wide (typically for base R plots) or long (typically for ggplot). With this, we can recreate the plot of transition probabilities shown above, in ggplot.

# Plot of transition probs as function of slope
plotData1 <- getPlotData(m = mod2, type = "tpm", format = "long")
ggplot(plotData1$slope, aes(slope, mle)) + 
    facet_wrap("prob") +
    geom_line() +
    geom_ribbon(aes(ymin = lci, ymax = uci), alpha = 0.3) +
    labs(y = "transition probability")

Similarly, we can create a plot of stationary state probabilities as functions of slope.

# Plot of stationary state probs as function of slope
plotData2 <- getPlotData(m = mod2, type = "stat", format = "long")
ggplot(plotData2$slope, aes(slope,  mle, col = factor(state))) + 
    geom_line() +
    geom_ribbon(aes(ymin = lci, ymax = uci, col = NULL,
                    fill = factor(state)), alpha = 0.3) +
    labs(fill = "state", col = "state", y = "stationary state probabilities")

Other features

new_data <- data.frame(temp = 0, slope = 17)
new_data

  temp slope
1    0    17

tpm <- predictTPM(m = mod2, newData = new_data, returnCI = TRUE)
tpm

$mle
, , 1

             toState1  toState2
fromState1 0.46504460 0.5349554
fromState2 0.03283051 0.9671695


$lci
, , 1

             toState1  toState2
fromState1 0.27371399 0.3327530
fromState2 0.01418787 0.9258727


$uci
, , 1

             toState1  toState2
fromState1 0.66724704 0.7262860
fromState2 0.07412731 0.9858121