Analyzing Masked Booby Movement Using HMM

Author

Affiliations

Marie-Pierre Etienne

ENSAI - CREST

https://marieetienne.github.io/MAF/

Published

December 9, 2024

library(tidyverse)

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.4     ✔ readr     2.1.5
✔ forcats   1.0.0     ✔ stringr   1.5.1
✔ ggplot2   3.5.1     ✔ tibble    3.2.1
✔ lubridate 1.9.3     ✔ tidyr     1.3.1
✔ purrr     1.0.2     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(ggpubr)
library(sf)

Linking to GEOS 3.12.1, GDAL 3.8.4, PROJ 9.4.0; sf_use_s2() is TRUE

library("rnaturalearth")
library(lubridate)
library(nlme)


Attaching package: 'nlme'

The following object is masked from 'package:dplyr':

    collapse

source(file = 'utils_HMM.R')
source(file = 'ICL.R')
library(moveHMM)

Loading required package: CircStats
Loading required package: MASS

Attaching package: 'MASS'

The following object is masked from 'package:dplyr':

    select

Loading required package: boot

library(depmixS4)

Loading required package: nnet
Loading required package: Rsolnp

Attaching package: 'depmixS4'

The following objects are masked from 'package:moveHMM':

    stationary, viterbi

library(circular)


Attaching package: 'circular'

The following objects are masked from 'package:CircStats':

    A1, A1inv, deg, I.0, I.1, I.p, pp.plot, rad, rose.diag, rstable

The following objects are masked from 'package:stats':

    sd, var

knitr::opts_chunk$set(echo = TRUE, comment = NA, cache = TRUE, message = FALSE,
                      warning = FALSE, eval = FALSE,
                      fig.align = "center")
options(ggplot2.discrete.colour=   scale_color_manual(values = wesanderson::wes_palette(name = "Darjeeling1")) )
couleur <-  wesanderson::wes_palette(name = "Darjeeling1")
scale2 <- function(x, na.rm = FALSE) (x - mean(x, na.rm = na.rm)) / ( sqrt((length(x)-1) / length(x)) *sd(x, na.rm) )

Data Presentation

We focus on studying the behavior of Red-footed Boobies (Sula sula).

Our dataset consists of three trajectories recorded for three Red-footed Boobies from Ilha do Meio, a Brazilian island part of the Fernando de Noronha archipelago in the Atlantic Ocean. 2017 04 23 SDAC05 MEIO BR004 trip 1 .csv, 2017 04 23 SDAC06 MEIO BR1705 trip 1 .csv, 2017 04 23 SDAC08 MEIO BR1703 trip 1 .csv,

fou_dta <- map_dfr(c('2017_04_23_SDAC05_MEIO_BR004_trip1.csv', 
              '2017_04_23_SDAC06_MEIO_BR1705_trip1.csv',
              '2017_04_23_SDAC08_MEIO_BR1703_trip1.csv'
              ),
                   function(f_){
                       ID <- f_ %>% 
                         stringr::str_split( ' ') %>% 
                       unlist() %>% stringr::str_subset('BR')
                     read_delim(file = file.path(f_ ),
                                delim = ";", escape_double = FALSE,
                                trim_ws = TRUE) %>%
                       dplyr::select(id, vol, datetime,  alt, lat, lon, dist.nid) %>%
                       rename(id_loc = id)  %>%
                       mutate(ID = ID) %>% 
                       rowid_to_column(var = "time_step")
                   })  %>% 
  mutate(dist_scaled = scale(dist.nid),
         dist_scaled_sq = dist_scaled^2,
         alt_scaled = scale(alt))

zone_box <- c(-32.75, -4.7, -31.8, -3.75)
zone_map <- ggmap::get_stamenmap(bbox = zone_box, 
                                 zoom = 10)
traj_plot <- ggmap::ggmap(zone_map) + 
    labs(x = "Longitude", y = "Latitude") +
    geom_path(data = fou_dta, aes(x = lon, y = lat, color = ID)) +
  theme(legend.position = "none", axis.text = element_blank(),
        axis.ticks = element_blank(),
        axis.title = element_blank())

world <- ne_countries(scale = "medium", returnclass = "sf") %>%
  st_crop(c(xmin = -90, ymin = -60, xmax = -30, ymax = 15))

zone_contour <-  fou_dta %>% 
  st_as_sf(coords = c("lon", "lat") ) %>% 
  st_set_crs(4326) %>% 
  st_bbox() %>% 
  st_as_sfc()

world_plot <- world  %>%  
  ggplot() + 
  geom_sf(fill = "#c9d0a3") + 
  geom_sf(data = zone_contour, fill = "white", color = "red") +
  theme(panel.background = element_rect(fill = "#99b3cc"))

gridExtra::grid.arrange(world_plot, traj_plot, nrow = 1)

Projection

The data for this booby are specified in latitude and longitude, which means they are based on angular coordinates on Earth’s surface. Methods described earlier rely on distance metrics (e.g., step lengths). While it’s possible to compute traveled distances using latitude and longitude, this involves formulas specific to spherical movements.

To circumvent this and use Euclidean distance, data are often projected into a local system. This projection must match the region of interest. Here, we use UTM coordinates for Zone 25 South. Geographic data manipulation is made easier in R using the sf package, which we use to project the data for easier distance calculations.

fou_dta_utm <- fou_dta %>% 
  st_as_sf(coords = c("lon", "lat")) %>% 
  mutate(dist_scaled = scale(dist.nid),
         dist_scaled_sq = dist_scaled^2,
         alt_scaled = scale(alt)) %>% 
  st_set_crs(4326)  %>% 
  st_transform(crs = 32725) %>% 
  mutate(Easting = st_coordinates(.)[,"X"],
         Northing = st_coordinates(.)[,"Y"])

Smoothing the Data

Data here are acquired at high frequency (one point every 10 seconds). Although GPS data have minimal error, the close temporal spacing of observations can add erraticity to trajectories. Such erraticity is particularly visible in movement metrics used to detect different activities.

To correct for these errors, we define a simple Gaussian linear hidden Markov model, as discussed earlier. The hidden dynamics are modeled as a random walk with an unknown innovation matrix, and the observation model assumes a Gaussian distribution centered on the true position with an unknown diagonal variance. The parameters are estimated using the EM algorithm in the MARSS package, and Kalman smoothing is used to estimate the true trajectory.

A smoothed trajectory preview is shown in the figure below, illustrating how smoothing reduces trajectory irregularities. However, note that this transformation involves working with modified data, whose relevance may vary in specific cases.

fou_dta_utm_smoothed %>% 
  dplyr::select(Easting, Easting_smoothed, Northing, Northing_smoothed) %>% 
  unite(col = "Raw", Easting, Northing, sep = "-") %>% 
  unite(col = "Smoothed", Easting_smoothed, Northing_smoothed, sep = "-") %>% 
  gather(key = "Trajectory", value = "Coordinates", factor_key = TRUE) %>% 
  separate(col = Coordinates, into = c("Longitude", "Latitude"), sep = "-", convert = TRUE) %>% 
  mutate(Trajectory = factor(Trajectory, labels = c("Raw", "Smoothed"))) %>% 
  ggplot() +
  aes(x = Longitude, y = Latitude, color = Trajectory) +
  geom_path() +
  coord_cartesian(ylim = c(955, 956) * 1e4, xlim = c(575, 580) * 1e3) +
  theme(axis.text = element_blank(), axis.ticks = element_blank(), axis.title = element_blank())

Fitting a Three-State Model

First, we fit a three-state model. The discussion about choosing the number of states is addressed later.

Defining Metrics

We aim to identify different behavioral types within the trajectory, such as foraging (associated with rapid direction changes) and goal-directed travel (characterized by straighter trajectories). Relevant metrics here are turning angle and step length, as these are biologically meaningful for these birds.

We fit two models differing by the considered movement metrics (and their corresponding emission distributions):

Step Length/Angle: Metrics from Morales et al. (2004), using gamma distribution for step length and von Mises distribution for angles.
Bivariate Velocity: Metrics from Gurarie, Andrews, and Laidre (2009), using two independent normal distributions.

Initializing the Algorithm

These models without covariates have 18 parameters (12 for emission distributions and 6 for transition probabilities). Optimization in such a parameter space is sensitive to initial values. We use a simple k-means classification to find plausible initial parameters for the regimes.

moveHMM_data <- fou_dta_utm_smoothed %>% 
  as.data.frame() %>%
  rowid_to_column(var = "id_point") %>% 
  dplyr::select(ID, id_point, time_step, Easting, Northing, Easting_smoothed, Northing_smoothed, 
                alt_scaled, dist_scaled, dist_scaled_sq) %>% 
  moveHMM::prepData(type = "UTM", coordNames = c("Easting_smoothed", "Northing_smoothed"))

depmix_data <- moveHMM_data %>% 
  as_tibble() %>% 
  mutate(v_p = step * cos(angle), v_r = step * sin(angle)) %>%
  replace_na(list(angle = 0, v_r = 0)) %>% 
  mutate(v_p = ifelse(is.na(v_p), step, v_p))
n_times <- depmix_data %>% 
  group_by(ID) %>% 
  summarise(n_times = n()) %>% 
  pull(n_times)

library(depmixS4)
source("utils_HMM.R")
n_states <- 3
set.seed(123)
initial_model <- depmixS4::depmix(list(v_p ~ 1, v_r ~ 1), data = depmix_data, 
                                  nstates = n_states,
                                  family = list(gaussian(), gaussian()), ntimes = n_times,
                                  respstart = get_init_depmix(depmix_data, nbStates = n_states),
                                  initdata = rep(1/n_states, n_states),
                                  transition = ~ 1)
depmix_fit <- depmixS4::fit(initial_model, verbose = FALSE, emcontrol = em.control(crit = "relative"))

rank_vector_depmix <-  posterior(depmix_fit) %>% 
  dplyr::select(state) %>% 
  mutate(step = depmix_data$step) %>% 
  group_by(state) %>% 
  summarise(mean_step = mean(step, na.rm = TRUE)) %>% 
  arrange(state) %>% 
  pull(mean_step) %>% 
  rank()
depmix_states <- posterior(depmix_fit) %>% 
  rename(Predicted_state = state) %>% 
  mutate(Predicted_state = rank_vector_depmix[Predicted_state] %>% 
           factor(levels = 1:n_states)) %>% 
  bind_cols(depmix_data, .) %>% 
  as_tibble() %>% 
  mutate(metric = "Bivariate Velocity")

par_init_moveHMM <- get_init_moveHMM(dta = moveHMM_data, nbStates = n_states)
set.seed(123)
moveHMM_fit <- moveHMM::fitHMM(data = moveHMM_data, nbStates = n_states, 
                               stepPar0 = par_init_moveHMM$stepPar0,
                               anglePar0 = par_init_moveHMM$anglePar0,
                               formula = ~ 1)

rank_vector_moveHMM <-  tibble(state = moveHMM::viterbi(moveHMM_fit)) %>% 
  mutate(step = depmix_data$step) %>% 
  group_by(state) %>% 
  summarise(mean_step = mean(step, na.rm = TRUE)) %>% 
  arrange(state) %>% 
  pull(mean_step) %>% 
  rank()
moveHMM_states <- moveHMM::stateProbs(moveHMM_fit) %>% 
  as_tibble() %>% 
  mutate(Predicted_state = rank_vector_moveHMM[moveHMM::viterbi(moveHMM_fit)] %>% 
           factor(levels = 1:n_states)) %>% 
  mutate(metric = "Step Length/Angle")
ICL(moveHMM_fit)

estimated_states <- moveHMM_states %>% 
  bind_rows(depmix_states)

Characterizing Hidden States

For both packages used, HMM parameters are estimated via maximum likelihood, and the most probable sequence of hidden states is retrieved using the Viterbi algorithm. Hidden states correspond to velocity and turning angle distributions. In trajectory terms, they describe segments between two positions (separated by 10 seconds).

estimated_states %>% 
  group_by(ID, metric) %>%
  mutate(Next_East = lead(x), Next_North = lead(y)) %>% 
  ggplot(aes(x = x, y = y)) +
  geom_segment(aes(xend = Next_East, yend = Next_North, color = Predicted_state)) +
  facet_wrap(~ metric) +
  labs(color = "Predicted State") +
  theme(axis.title = element_blank(), axis.text = element_blank(), axis.ticks = element_blank())

References

Gurarie, Eliezer, Russel D Andrews, and Kristin L Laidre. 2009. “A Novel Method for Identifying Behavioural Changes in Animal Movement Data.” Ecology Letters 12 (5): 395–408.

Morales, Juan Manuel, Daniel T Haydon, Jacqui Frair, Kent E Holsinger, and John M Fryxell. 2004. “Extracting More Out of Relocation Data: Building Movement Models as Mixtures of Random Walks.” Ecology 85 (9): 2436–45.