Purpose of this workshop: We will use the seals’ dataset from Tutorial 2 to build an hidden Markov model

the R packages “moveHMM” and “momentuHMM” are the two main libraries we will use.
Remember, HMMs run on regular time series data! (our data from Tutorial 2 is already regular)

Resources:
McClintock, B. T., & Michelot, T. (2018). momentuHMM: R package for generalized hidden Markov models of animal movement. Methods in Ecology and Evolution, 9(6), 1518–1530. https://doi.org/10.1111/2041-210X.12995

Michelot, T., Langrock, R., & Patterson, T. A. (2016). moveHMM: An R package for the statistical modelling of animal movement data using hidden Markov models. Methods in Ecology and Evolution, 7(11), 1308–1315. https://doi.org/10.1111/2041-210X.12578

See also the vignettes:
https://cran.r-project.org/web/packages/momentuHMM/index.html
https://cran.r-project.org/web/packages/moveHMM/index.html

It will be up to you to install them and see if you have all the packages ready.
See that we clean the environment first and then start loading packages.

#clean environment
rm(list = ls())
gc()
##          used (Mb) gc trigger (Mb) max used (Mb)
## Ncells 523291 28.0    1163211 62.2   663811 35.5
## Vcells 963071  7.4    8388608 64.0  1817995 13.9
library(readr)
library(moveHMM)
library(momentuHMM)

Load the data and:
1. make sure you have the name of the columns as the package requires, here we need a column called “ID” - capital letters
2. we run the prepData function

SealsData<- read_csv("~/MTP_Practicals/Seals_GPS_Data_MTP.csv")
SealsData<-as.data.frame(SealsData)#we turn the dataset in a simple dataset format
head(SealsData)
##                  date    id mean_dive_depth diving_duration_sec       x       y
## 1 2017-03-06 21:36:07 dewey               0                   0 1147644 6474366
## 2 2017-03-06 21:56:07 dewey               0                   0 1147644 6474363
## 3 2017-03-06 22:16:07 dewey               0                   0 1147643 6474361
## 4 2017-03-06 22:36:07 dewey               0                   0 1147643 6474358
## 5 2017-03-06 22:56:07 dewey               0                   0 1147642 6474356
## 6 2017-03-06 23:16:07 dewey               0                   0 1147642 6474353
colnames(SealsData)[2]<-"ID"  #the package needs the column to be called "ID"
SealsData$ID<-as.factor(SealsData$ID) 

#see what prepData does from the R help and inspect the results and the tracks.
data_prepped <- prepData(SealsData, type = 'UTM', coordNames = c("x", "y"))
head(data_prepped)
##      ID     step         angle                date mean_dive_depth
## 1 dewey 2.653222            NA 2017-03-06 21:36:07               0
## 2 dewey 2.653222 -4.866972e-07 2017-03-06 21:56:07               0
## 3 dewey 2.653222 -4.865398e-07 2017-03-06 22:16:07               0
## 4 dewey 2.653222 -4.821456e-07 2017-03-06 22:36:07               0
## 5 dewey 2.653222 -4.867684e-07 2017-03-06 22:56:07               0
## 6 dewey 3.829366  2.798740e+00 2017-03-06 23:16:07               0
##   diving_duration_sec       x       y
## 1                   0 1147644 6474366
## 2                   0 1147644 6474363
## 3                   0 1147643 6474361
## 4                   0 1147643 6474358
## 5                   0 1147642 6474356
## 6                   0 1147642 6474353
#let's discuss what you see! also, go slowly, you will hit the enter button for each plot
plot(data_prepped)

Set the model and select starting parameters.

#we go from m to Km - helps the numerical stability of the model
data_prepped$step<-data_prepped$step/1000

#we ask for 2 behavioural states
nbStates <- 2 
stateNames <- c("state1", "state2")

#which distributions?
dist <- list(step = "gamma", angle = "vm")

#which parameters?
stepPar0 <- c(0.2, 1, 0.01, 0.5) 
anglePar0 <- c(0, 0, 1,1) # angle mean and concentration parameters for each state
#high concentration represents consistent direction, low concentration parameter indicates more erratic changes in direction

Run the model and display results

fitSeals <- fitHMM(data = data_prepped, nbStates = 2,dist=dist,
                   Par0=list(step=stepPar0,angle=anglePar0), 
                   formula = ~1,
                   estAngleMean = list(angle=TRUE))

#plot results, which ecological interpretation can you give to state 1 and state 2?
plot(fitSeals)
## Decoding state sequence... DONE

#see model summary, what is transition probability matrix telling you?
fitSeals
## Value of the maximum log-likelihood: -2290.295 
## 
## 
## step parameters:
## ----------------
##        state 1   state 2
## mean 0.1304798 0.8087957
## sd   0.1467984 0.4121233
## 
## angle parameters:
## -----------------
##                   state 1     state 2
## mean          -0.02641102 -0.01123537
## concentration  0.65133783  1.88135018
## 
## Regression coeffs for the transition probabilities:
## ---------------------------------------------------
##                1 -> 2    2 -> 1
## (Intercept) -2.607925 -2.000379
## 
## Transition probability matrix:
## ------------------------------
##           state 1   state 2
## state 1 0.9313699 0.0686301
## state 2 0.1191631 0.8808369
## 
## Initial distribution:
## ---------------------
##      state 1      state 2 
## 0.9999881833 0.0000118167

The viterbi algorithm is used to estimate the behavioural states, and we attach them to our dataset

data_prepped$states <- viterbi(fitSeals)

#we also see how many datapoints belong to each state, are the animals spending more time in state 1 or state 2?
table(data_prepped$states)
## 
##    1    2 
## 1140  657

We check our model

# compute the pseudo-residuals, here you will see the Markov chain, and the residual autocorrelation.
# is there residual autocorrelation left?
pr <- pseudoRes(fitSeals) 
# time series, qq-plots, and ACF of the pseudo-residuals 
plotPR(fitSeals)
## Computing pseudo-residuals...

# here we visualize the probability of each data point belonging to a specific behavioral state.
# Discuss with us what you see in these plots.
# try with other animal ids
plotStates(fitSeals,animals="dewey")
## Decoding state sequence... DONE
## Computing state probabilities... DONE