Capuchins Movement Segmentation

AKDE Index Daily Travels

We applied Hidden Markov Models (HMMs; McClintock et al., 2020) to analyze the movement patterns of black capuchins in plantation landscapes and in INP. This approach allowed us to identify and compare discrete movement states between these contrasting landscapes and between contrasting environments within plantation landscapes (pine stands vs forest fragments). Within plantation landscapes we evaluated if movement states showed seasonal variation between native forest patches and pine stands associated to the seasonal bark-stripping behavior of capuchins (Mikich and Liebsch, 2014), expecting shorter steps and more convoluted movement patterns associated with this activity. The dataset consisted of GPS locations recorded at 30-minute intervals throughout the activity period (6:30 AM–8:30 PM) for four groups in plantations and four groups in INP. Data spanned multiple years (2008, 2009, 2010, 2022, and 2023) covering the months of March to December, with June excluded due to inconsistent data collection. We classified movement, based on step length (distance between consecutive GPS locations) and turning angle (change in direction), into two behavioral states: Resident (short distances, sharp turns) and Transit (long distances, small turns). Step length represents movement magnitude (speed), while turning angle indicates directional persistence (Langrock et al., 2012). Resident behavior, linked to resting, social interactions, or slow foraging, had an initial mean step length of 20 m (SD = 20 m) and turning angle concentration (k) of 0.1. Transit behavior, associated with direct travel between resources, had a mean step length of 300 m (SD = 300 m) and k = 0.5. We employed the Viterbi algorithm to decode the most probable sequence of behavioral states for each individual (Joo et al., 2013). We evaluated model convergence by verifying stable parameter estimates across multiple runs with varying initial values, ensuring that the global maximum of the likelihood function was identified (McClintock and Michelot, 2018). These analyses were conducted using the R package ‘momentuHMM’ (McClintock and Michelot, 2018).

1 R Packages

library(sf)
library(sp)
library(amt)
library(tidyverse)
library(momentuHMM)
library(ggplot2)
library(conflicted)
library(patchwork)
library(flextable)
library(gridExtra)
library(grid)
library(showtext)
font_add("Calibri", "C:/Windows/Fonts/calibri.ttf")  
font_add("Calibri Bold", "C:/Windows/Fonts/calibrib.ttf")  
font_add("Calibri Italic", "C:/Windows/Fonts/calibrii.ttf")
conflict_prefer("filter", "dplyr")
conflict_prefer("select", "dplyr") 
showtext_auto()

1.1 Plantations processed data

hmm <- bind_rows(lapply(c('tilo_hmm.RDS', 'toxo_hmm.RDS', 'jasy_hmm.RDS', 'col_palomo.RDS'), readRDS))

Extracting the information from Timber company Shape file

shp <- 
  st_read("UNION.shp") %>%
  mutate(tipouso_grouped = ifelse(
    tipouso == "BNSM",
    "BNSM",
    "plantations")) %>% 
  mutate(sensibilidad = case_when(
    tipouso == "PITA" & apl >= 2012 & apl <= 2017 ~ "sensible",
    TRUE ~ "insensible"))

crs_shp <- st_crs(shp)

shp2 <- 
  st_read("UNION.shp") %>%
  mutate(tipouso_grouped = ifelse(
    tipouso == "BNSM",
    "BNSM",
    "plantations")) %>% 
  mutate(ages =case_when(
      tipouso_grouped == "plantations" & apl >= 2018 ~ "young",
      tipouso_grouped == "plantations" & apl >= 2013 & apl <= 2017 ~ "medium",
      tipouso_grouped == "plantations" & apl <= 2012 ~ "old",
      TRUE ~ NA_character_  
    ))


shp <- shp %>%
  mutate(tipouso_grouped = ifelse(tipouso == "BNSM", "BNSM", "plantations"))

These will allow me to calculate minimum distantes to different landscape features

shp_bnsm <- shp %>% filter(tipouso == "BNSM")
shp_medium <- shp2 %>% filter(ages == 'medium')
shp_old<- shp2 %>% filter(ages == 'old')
shp_young <- shp2 %>% filter(ages == 'young')

hmm <- 
 bind_rows(lapply(c(
   'tilo_hmm.RDS',
   'toxo_hmm.RDS',
   'jasy_hmm.RDS',
   'col_palomo.RDS'),
   readRDS)) %>%
  rename(lon = x_, lat = y_, t=t_) %>% 
  select(lat, lon, t, id)


dist_sensible <- hmm_sf2[,"dist_sensible"] %>% 
  st_drop_geometry() %>% 
  pull(dist_sensible) %>% 
  as.numeric()

dist_bnsm <- hmm_sf2[,"dist_bnsm"] %>% 
  st_drop_geometry() %>% 
  pull(dist_bnsm) %>% 
  as.numeric()

dist_old <- hmm_sf2[,"dist_old"] %>% 
  st_drop_geometry() %>% 
  pull(dist_old) %>% 
  as.numeric()

dist_young <- hmm_sf2[,"dist_young"] %>% 
  st_drop_geometry() %>% 
  pull(dist_young) %>% 
  as.numeric()


hmm <- hmm %>% 
   mutate(
    X = X,  
    Y = Y,  
    dist_bnsm = dist_bnsm,  
    dist_sensible = dist_sensible,
    dist_young = dist_young,
    dist_old = dist_old)

hmm_sf <- hmm %>% 
  st_as_sf(coords = c('lon', 'lat'),
           crs = 4326)



if (st_crs(hmm_sf) != st_crs(shp)) {
  shp <- st_transform(shp, st_crs(hmm_sf))  
}


hmm_sf_joined <- st_join(hmm_sf, shp[, c("tipouso", "apl")])

hmm_sf_joined <- hmm_sf_joined %>%
  mutate(tipouso_grouped = ifelse(tipouso == "BNSM", "BNSM", "plantations"))

hmm_sf_joined <- hmm_sf_joined %>%
  mutate(
    tipouso_grouped = case_when(
      tipouso %in% c("BNSM", "FCAM") ~ "BNSM",  
      TRUE ~ "plantations"                    
    ),
    
    apl_grouped = case_when(
      tipouso_grouped == "plantations" & apl >= 2017 ~ "young",
      tipouso_grouped == "plantations" & apl >= 2013 & apl <= 2016 ~ "medium",
      tipouso_grouped == "plantations" & apl <= 2012 ~ "old",
      TRUE ~ NA_character_  
    ))

hmm <- hmm_sf_joined %>%
  st_drop_geometry() %>%  
  bind_cols(st_coordinates(hmm_sf_joined) %>% 
              as.data.frame() %>% 
              rename(lon = X, lat = Y))

Explorative plots

hmm_sf <- hmm %>% 
  st_as_sf(coords = c('lon', 'lat'),
           crs = 4326)
datos_filtrados <- hmm_sf %>%
  mutate(
    dist_sensible = as.numeric(dist_sensible),  
    dist_bnsm = as.numeric(dist_bnsm)) %>%
  pivot_longer(
    cols = c(dist_sensible, dist_bnsm),
    names_to = "metric",
    values_to = "distancia")

Constrasting distances to plantations and native forest fragments

p_col<-
  datos_filtrados %>% filter (id == 'col', metric == 'dist_sensible') %>% 
  ggplot() +
  geom_sf(aes(color = distancia),
          size = 3) +
  scale_color_gradientn(
    colors = c("gold", "darkorange", "orangered", "red1", "firebrick3", "darkred"),
    name = "Distance (m)") +
  labs(
    x = "Longitude",
    y = "Latitud") +
  theme_minimal() +
  theme(
    panel.background = element_rect(fill = "white", color = NA),  
    plot.background = element_rect(fill = "white", color = NA),
    legend.position = "none",
    text = element_text(family = "Calibri", size = 34),
    axis.text.x = element_text(angle = 45, hjust = 1),
    plot.margin = margin(0, 0, 0, 0))

p_tilo<-
  datos_filtrados %>% filter (id == 'tilo', metric == 'dist_sensible') %>% 
  ggplot() +
  geom_sf(aes(color = distancia),
          size = 3) +
  scale_color_gradientn(
    colors = c("gold", "darkorange", "orangered", "red1", "firebrick3", "darkred"),
    name = "Distance (m)") +
  labs(
    x = "Longitude",
    y = " ") +
  theme_minimal() +
  theme(
    panel.background = element_rect(fill = "white", color = NA),  
    plot.background = element_rect(fill = "white", color = NA),
    legend.position = "right",
    text = element_text(family = "Calibri", size = 34),
    axis.text.x = element_text(angle = 45, hjust = 1),
    plot.margin = margin(0, 0, 0, 0))

p_jasy<-
  datos_filtrados %>% filter (id == 'jasy', metric == 'dist_sensible') %>% 
  ggplot() +
  geom_sf(aes(color = distancia),
          size = 3) +
  scale_color_gradientn(
    colors = c("gold", "darkorange", "orangered", "red1", "firebrick3", "darkred"),
    name = "Distance (m)") +
  labs(
    x = "Longitude",
    y = " ") +
  theme_minimal() +
  theme(
    panel.background = element_rect(fill = "white", color = NA),  
    plot.background = element_rect(fill = "white", color = NA),
    legend.position = "none",
    text = element_text(family = "Calibri", size = 34),
    axis.text.x = element_text(angle = 45, hjust = 1),
    plot.margin = margin(0, 0, 0, 0))

p_toxo<-
  datos_filtrados %>% filter (id == 'toxo', metric == 'dist_sensible') %>% 
  ggplot() +
  geom_sf(aes(color = distancia),
          size = 3) +
  scale_color_gradientn(
    colors = c("gold", "darkorange", "orangered", "red1", "firebrick3", "darkred"),
    name = "Distance (m)") +
  labs(
    x = "Longitude",
    y = "Latitud") +
  theme_minimal() +
  theme(
    panel.background = element_rect(fill = "white", color = NA),  
    plot.background = element_rect(fill = "white", color = NA),
    legend.position = "none",
    text = element_text(family = "Calibri", size = 34),
    axis.text.x = element_text(angle = 45, hjust = 1),
    plot.margin = margin(0, 0, 0, 0))

Native Forest Fragments

p_col_nf<-
  datos_filtrados %>% filter (id == 'col', metric == 'dist_bnsm') %>% 
  ggplot() +
  geom_sf(aes(color = distancia),
          size = 3) +
  scale_color_gradientn(
    colors = c("darkolivegreen","darkolivegreen3", "darkolivegreen2", "darkolivegreen1", "palegreen1"),
    name = "Distance (m)"
  )  +
  labs(
    x = "Longitude",
    y = "Latitud") +
  theme_minimal() +
  theme(
    panel.background = element_rect(fill = "white", color = NA),  
    plot.background = element_rect(fill = "white", color = NA),
    legend.position = "none",
    text = element_text(family = "Calibri", size = 34),
    axis.text.x = element_text(angle = 45, hjust = 1),
    plot.margin = margin(0, 0, 0, 0))

p_tilo_nf<-
  datos_filtrados %>% filter (id == 'tilo', metric == 'dist_bnsm') %>% 
  ggplot() +
  geom_sf(aes(color = distancia),
          size =3) +
  scale_color_gradientn(
    colors = c("darkolivegreen","darkolivegreen3", "darkolivegreen2", "darkolivegreen1", "palegreen1"),
    name = "Distance (m)"
  )  +
  labs(
    x = "Longitude",
    y = " ") +
  theme_minimal() +
  theme(
    panel.background = element_rect(fill = "white", color = NA),  
    plot.background = element_rect(fill = "white", color = NA),
    legend.position = "right",
    text = element_text(family = "Calibri", size = 34),
    axis.text.x = element_text(angle = 45, hjust = 1),
    plot.margin = margin(0, 0, 0, 0))

p_jasy_nf<-
  datos_filtrados %>% filter (id == 'jasy', metric == 'dist_bnsm') %>% 
  ggplot() +
  geom_sf(aes(color = distancia),
          size = 3) +
  scale_color_gradientn(
    colors = c("darkolivegreen","darkolivegreen3", "darkolivegreen2", "darkolivegreen1", "palegreen1"),
    name = "Distance (m)"
  )  +
  labs(
    x = "Longitude",
    y = " ") +
  theme_minimal() +
  theme(
    panel.background = element_rect(fill = "white", color = NA),  
    plot.background = element_rect(fill = "white", color = NA),
    legend.position = "none",
    text = element_text(family = "Calibri", size = 34),
    axis.text.x = element_text(angle = 45, hjust = 1),
    plot.margin = margin(0, 3, 0, 0))

p_toxo_nf<-
  datos_filtrados %>% filter (id == 'toxo', metric == 'dist_bnsm') %>% 
  ggplot() +
  geom_sf(aes(color = distancia),
          size = 3) +
  scale_color_gradientn(
    colors = c("darkolivegreen","darkolivegreen3", "darkolivegreen2", "darkolivegreen1", "palegreen1"),
    name = "Distance(m)"
  )  +
  labs(
    x = "Longitude",
    y = "Latitud") +
  theme_minimal() +
  theme(
    panel.background = element_rect(fill = "white", color = NA),  
    plot.background = element_rect(fill = "white", color = NA),
    legend.position = "none",
    text = element_text(family = "Calibri", size = 34),
    axis.text.x = element_text(angle = 45, hjust = 1),
    plot.margin = margin(0, 0, 0, 0))

p1 <-  ((p_col_nf | p_tilo_nf) / (p_toxo_nf | p_jasy_nf)) +
  plot_layout(guides = "collect") & 
  theme(plot.background = element_rect(fill = "white", color = NA))

p2 <- ((p_col | p_tilo) / (p_toxo | p_jasy)) +
  plot_layout(guides = "collect") & 
  theme(plot.background = element_rect(fill = "white", color = NA))

p1 + plot_annotation(title = "Distances to Native Forest")& theme(
  plot.title = element_text(size = 40, face = "bold"))

p2 +plot_annotation(
  title = "Distance to Medium-aged Pine Stands") & theme(
  plot.title = element_text(size = 40, face = "bold"))

2 Preparing data for HMM

There are some zeros, so I continue with the ZeroMass included.

I had to remove NAs, likely due to some issues when distances were calculated from the shape file. It was about 30 points, nothing serious. I also scaled both distances because many Warnings were triggered due to the unusual distribution of distances… many zeros and some values of 1500 or more. I keep the columns with the variables in their original scale.

hmm_prep <- 
  prepData(
    hmm, type = "UTM",
    coordNames = c("x","y"))

cutoff_date <- as.POSIXct("2022-07-21", tz = "UTC")  # Ajusta el año según corresponda

low <- hmm_prep %>% filter(t < cutoff_date)
high <- hmm_prep %>% filter(t >= cutoff_date)

3 Initial parameters

stepMean0_2 <- c(m1 = 20, m2 = 300)
stepSD0_2 <- c(sd1 = 20, sd2 = 300)
zeromass_2 <- c(zm1 = 0.05, zm2 =0.05) 
angleCon0_2 <- c(rho1  = 0.1, rho2 = 0.5)

stateNames_2 <- c("resident","transit")

dist_2 <- list(step = "gamma", angle = "wrpcauchy")
Par0_2 <- list(step=c(stepMean0_2, stepSD0_2), angle = c(angleCon0_2))

4 Model Fitting

4.1 Full year Models

Full-year models include data from March to December, that is low and high bark-stripping seasons. ‘Sens’ stands for medium-aged pine stands

# Model 1: Null model
hmm2_null <-
  fitHMM(
    data = hmm_prep,
    nbStates = 2, dist = dist_2,
    Par0 = Par0_2, stateNames = stateNames_2,
    formula = ~ 1)

# Model 2: Medium-aged distance
hmm2_sens <-
  fitHMM(
    data = hmm_prep,
    nbStates = 2, dist = dist_2,
    Par0 = Par0_2, stateNames = stateNames_2,
    formula = ~ dist_sensible)

# Model 4: Old distance
hmm2_old <-
  fitHMM(
    data = hmm_prep,
    nbStates = 2, dist = dist_2,
    Par0 = Par0_2, stateNames = stateNames_2,
    formula = ~ dist_old)

# Model 5: Young distance
hmm2_young <-
  fitHMM(
    data = hmm_prep,
    nbStates = 2, dist = dist_2,
    Par0 = Par0_2, stateNames = stateNames_2,
    formula = ~ dist_young)

# Model 6: Medium-aged + BN+SM distance
hmm2_sens_bnsm <-
  fitHMM(
    data = hmm_prep,
    nbStates = 2, dist = dist_2,
    Par0 = Par0_2, stateNames = stateNames_2,
    formula = ~ dist_sensible + dist_bnsm)

# Model 7: Medium-aged + Old distance
hmm2_old_sens <-
  fitHMM(
    data = hmm_prep,
    nbStates = 2, dist = dist_2,
    Par0 = Par0_2, stateNames = stateNames_2,
    formula = ~ dist_sensible + dist_old)

# Model 8: Medium-aged + Young distance
hmm2_young_sens <-
  fitHMM(
    data = hmm_prep,
    nbStates = 2, dist = dist_2,
    Par0 = Par0_2, stateNames = stateNames_2,
    formula = ~ dist_sensible + dist_young)

# Model 9: Full model
hmm2_full <-
  fitHMM(
    data = hmm_prep,
    nbStates = 2, dist = dist_2,
    Par0 = Par0_2, stateNames = stateNames_2,
    formula = ~ dist_sensible + dist_bnsm + dist_old + dist_young)

Failed fitting of the model including Distance to young and medium-aged stands, optimization algorithm nlm supplied a non-finite value. Model specification leads to numerical instability, possibly due to collinearity between covariates. Consequently, this model is not feasible to run as currently defined, and thus excluded.

4.1.1 Full year AIC comparison

aic_year <- AIC(hmm2_null, hmm2_sens, hmm2_bnsm, hmm2_old, hmm2_young, 
                hmm2_sens_bnsm, hmm2_old_sens, hmm2_full)

aic_year_df <- as.data.frame(aic_year) %>%
  mutate(Model = c("Forest Fragment", "Young", "Medium", "Null", 
                   "Old", "Old Pine + Medium", "Medium + Forest Fragment", "Full Model"),
         deltaAIC = AIC - min(AIC)) %>%
  arrange(AIC) %>%
  select(Model, AIC, deltaAIC)

aic_year_table <- flextable(aic_year_df) %>%
  colformat_num(j = c("AIC", "deltaAIC"), digits = 3) %>%
  theme_tron_legacy() %>%
  autofit()

aic_year_table

Model	AIC	deltaAIC
Forest Fragment	81,785.64	0.0000
Young	81,935.88	150.2406
Medium	81,936.67	151.0253
Null	81,937.87	152.2320
Old	81,938.75	153.1074
Old Pine + Medium	82,294.64	508.9998
Medium + Forest Fragment	83,084.88	1,299.2429
Full Model	83,096.46	1,310.8182

hmm_2states <- viterbi(hmm2_bnsm)
prop.table(table(hmm_2states))

## hmm_2states
##       1       2 
## 0.74465 0.25535

hmm2_bnsm

## Value of the maximum log-likelihood: -40881.82 
## 
## 
## step parameters:
## ----------------
##      resident  transit
## mean 76.10556 229.8218
## sd   55.17141 148.0904
## 
## angle parameters:
## -----------------
##                resident   transit
## mean          0.0000000 0.0000000
## concentration 0.1464638 0.5154257
## 
## Regression coeffs for the transition probabilities:
## ---------------------------------------------------
##                  1 -> 2      2 -> 1
## (Intercept) -2.23279017 -0.47837146
## dist_bnsm    0.01445949 -0.01344355
## 
## Transition probability matrix (based on mean covariate values):
## ---------------------------------------------------------------
##           resident   transit
## resident 0.8687853 0.1312147
## transit  0.3107032 0.6892968
## 
## Initial distribution:
## ---------------------
##  resident   transit 
## 0.8068209 0.1931791

4.2 Low Bark-stripping models

Only fitting HMM for data collected during low bark-stripping month, see Figure B.2 panel A for more details

hmm2_low_null <-
  fitHMM(
    data = low,
    nbStates = 2, dist = dist_2,
    Par0 = Par0_2, stateNames = stateNames_2,
    formula = ~ 1)

hmm2_low_full <-
  fitHMM(
    data = low,
    nbStates = 2, dist = dist_2,
    Par0 = Par0_2, stateNames = stateNames_2,
    formula = ~ dist_bnsm + dist_sensible + dist_young + dist_old)

hmm2_low_sens <-
  fitHMM(
    data = low,
    nbStates = 2, dist = dist_2,
    Par0 = Par0_2, stateNames = stateNames_2,
    formula = ~ dist_sensible)

hmm2_low_bnsm <-
  fitHMM(
    data = low,
    nbStates = 2, dist = dist_2,
    Par0 = Par0_2, stateNames = stateNames_2,
    formula = ~ dist_bnsm)

hmm2_low_sens_bnsm <-
  fitHMM(
    data = low,
    nbStates = 2, dist = dist_2,
    Par0 = Par0_2, stateNames = stateNames_2,
    formula = ~ dist_sensible+ dist_bnsm)

hmm2_low_young <-
  fitHMM(
    data = low,
    nbStates = 2, dist = dist_2,
    Par0 = Par0_2, stateNames = stateNames_2,
    formula = ~ dist_young)

hmm2_low_old<-
  fitHMM(
    data = low,
    nbStates = 2, dist = dist_2,
    Par0 = Par0_2, stateNames = stateNames_2,
    formula = ~ dist_old)

hmm2_low_young_sens_bnsm <-
  fitHMM(
    data = low,
    nbStates = 2, dist = dist_2,
    Par0 = Par0_2, stateNames = stateNames_2,
    formula = ~ dist_sensible + dist_young + dist_bnsm)

hmm2_low_old_sens_bnsm<-
  fitHMM(
    data = low,
    nbStates = 2, dist = dist_2,
    Par0 = Par0_2, stateNames = stateNames_2,
    formula = ~ dist_sensible + dist_old + dist_bnsm)

Failed fitting of the model including Distance to native forest, old and medium-aged stands, optimization algorithm nlm supplied a non-finite value. Model specification leads to numerical instability, possibly due to collinearity between covariates. Consequently, this model is not feasible to run as currently defined, and thus excluded.

4.2.1 Low bark-stripping AIC comparison

aic_low <-AIC(hmm2_low_null, hmm2_low_sens, hmm2_low_sens_bnsm, hmm2_low_old ,hmm2_low_young, hmm2_low_young_sens_bnsm, hmm2_low_old_sens_bnsm)

aic_low_df <- as.data.frame(aic_low) %>%
  mutate(Model = c("Medium + Forest Fragment", "Medium + Young + Forest Fragment", "Medium", "Null", 
                   "Young","Old", "Old + Medium + Forest Fragment"),
         deltaAIC = AIC - min(AIC)) %>%
  arrange(AIC) %>%
  select(Model, AIC, deltaAIC)

aic_low_table <- flextable(aic_low_df) %>%
  colformat_num(j = c("AIC", "deltaAIC"), digits = 3) %>%
  theme_tron_legacy() %>%
  autofit()

aic_low_table

Model	AIC	deltaAIC
Medium + Forest Fragment	28,484.91	0.000000
Medium + Young + Forest Fragment	28,489.30	4.387503
Medium	28,492.08	7.171087
Null	28,512.33	27.425731
Young	28,512.34	27.434699
Old	28,513.07	28.163874
Old + Medium + Forest Fragment	28,514.87	29.966900

hmm_2states_low <- viterbi(hmm2_low_sens_bnsm)
prop.table(table(hmm_2states_low))

## hmm_2states_low
##         1         2 
## 0.5347187 0.4652813

hmm2_low_sens_bnsm

## Value of the maximum log-likelihood: -14231.65 
## 
## 
## step parameters:
## ----------------
##      resident   transit
## mean 50.08417 148.76993
## sd   33.12291  98.58646
## 
## angle parameters:
## -----------------
##                resident   transit
## mean          0.0000000 0.0000000
## concentration 0.1063557 0.4054252
## 
## Regression coeffs for the transition probabilities:
## ---------------------------------------------------
##                     1 -> 2        2 -> 1
## (Intercept)   -1.309465231 -0.9496168657
## dist_sensible -0.000700878 -0.0003953991
## dist_bnsm      0.022186113 -0.0370060607
## 
## Transition probability matrix (based on mean covariate values):
## ---------------------------------------------------------------
##           resident   transit
## resident 0.7718309 0.2281691
## transit  0.1805468 0.8194532
## 
## Initial distribution:
## ---------------------
##  resident   transit 
## 0.4583366 0.5416634

pr_low <- pseudoRes(hmm2_low_sens_bnsm)
qqnorm(pr_low$stepRes,
       pch = 21,          
       col = "midnightblue",   
       cex = 2.2,         
       main = "Q-Q Plot of Step Residuals",
       xlab = "Theoretical Quantiles",
       ylab = "Sample Quantiles",
       cex.main = 2.2,      
       cex.lab = 1.8,     
       cex.axis = 1.6)

qqnorm(pr_low$angleRes,
       pch = 21,          
       col = "midnightblue",   
       cex = 2.2,         
       main = "Q-Q Plot of Turning angles Residuals",
       xlab = "Theoretical Quantiles",
       ylab = "Sample Quantiles",
       cex.main = 2.2,      
       cex.lab = 1.8,     
       cex.axis = 1.6)

4.2.2 Stationary Plots (Low)

plotStationary(
  hmm2_low_sens_bnsm, 
  plotCI = TRUE, 
  col = c("#0072B2", "#d33420"),  
  cex.axis = 2, 
  cex.lab = 2,
  cex.legend = 1.8,
  cex.main = 2,
  lwd =2,
  legend.pos = "topright")

Figure 3. Predicted stationary probabilities for the movement states (Resident and Transit) of black capuchin monkeys (Sapajus nigritus) in pine plantations, derived from two-state Hidden Markov Models. Each line represents the mean probability of a state, and the shaded ribbons indicate 95% confidence intervals. The probabilities are shown as a function of the minimum distance to young, medium-aged pine stands and native forest fragments, during low bark-stripping months (March to May).

4.3 High Bark-stripping Models

hmm2_high_null <-
  fitHMM(
    data = high,
    nbStates = 2, dist = dist_2,
    Par0 = Par0_2, stateNames = stateNames_2,
    formula = ~ 1)

hmm2_high_sens_bnsm <-
  fitHMM(
    data = high,
    nbStates = 2, dist = dist_2,
    Par0 = Par0_2, stateNames = stateNames_2,
    formula = ~ dist_sensible + dist_bnsm)

hmm2_high_sens <-
  fitHMM(
    data = high,
    nbStates = 2, dist = dist_2,
    Par0 = Par0_2, stateNames = stateNames_2,
    formula = ~ dist_sensible)

hmm2_high_full <-
  fitHMM(
    data = high,
    nbStates = 2, dist = dist_2,
    Par0 = Par0_2, stateNames = stateNames_2,
    formula = ~ dist_sensible + dist_bnsm + dist_old + dist_young)

hmm2_high_young <-
  fitHMM(
    data = high,
    nbStates = 2, dist = dist_2,
    Par0 = Par0_2, stateNames = stateNames_2,
    formula = ~ dist_young)

hmm2_high_old<-
  fitHMM(
    data = high,
    nbStates = 2, dist = dist_2,
    Par0 = Par0_2, stateNames = stateNames_2,
    formula = ~ dist_old)

hmm2_high_young_sens_bnsm <-
  fitHMM(
    data = high,
    nbStates = 2, dist = dist_2,
    Par0 = Par0_2, stateNames = stateNames_2,
    formula = ~ dist_sensible + dist_bnsm + dist_young)

hmm2_high_old_sens_bnsm<-
  fitHMM(
    data = high,
    nbStates = 2, dist = dist_2,
    Par0 = Par0_2, stateNames = stateNames_2,
    formula = ~ dist_sensible + dist_bnsm + dist_old)

4.3.1 High bark-stripping AIC comparison

aic_high <-AIC(hmm2_high_null, hmm2_high_sens_bnsm, hmm2_high_sens, hmm2_high_old, hmm2_high_young,
               hmm2_high_young_sens_bnsm, hmm2_high_old_sens_bnsm)

aic_high_df <- as.data.frame(aic_high) %>%
  mutate(Model = c("Medium + Young + Forest Fragment", "Medium + Old + Forest Fragment", "Medium + Forest Fragment", "Medium", "Null","Young", "Old"),
         deltaAIC = AIC - min(AIC)) %>%
  arrange(AIC) %>%
  select(Model, AIC, deltaAIC)

aic_high_table <- flextable(aic_high_df) %>%
  colformat_num(j = c("AIC", "deltaAIC"), digits = 3) %>%
  theme_tron_legacy() %>%
  autofit()

aic_high_table

Model	AIC	deltaAIC
Medium + Young + Forest Fragment	53,133.34	0.00000
Medium + Old + Forest Fragment	53,150.00	16.66230
Medium + Forest Fragment	53,151.07	17.73671
Medium	53,260.34	127.00547
Null	53,266.35	133.01297
Young	53,267.04	133.70222
Old	53,271.10	137.76008

hmm_2states_high <- viterbi(hmm2_high_young_sens_bnsm)

hmm2_high_young_sens_bnsm

## Value of the maximum log-likelihood: -26551.67 
## 
## 
## step parameters:
## ----------------
##      resident  transit
## mean 87.30877 281.0000
## sd   62.79430 169.3918
## 
## angle parameters:
## -----------------
##                resident   transit
## mean          0.0000000 0.0000000
## concentration 0.1625826 0.5663035
## 
## Regression coeffs for the transition probabilities:
## ---------------------------------------------------
##                      1 -> 2        2 -> 1
## (Intercept)   -2.5168271847 -0.3306755996
## dist_sensible  0.0008183783 -0.0007597368
## dist_bnsm      0.0090843200 -0.0194060404
## dist_young    -0.0002480793  0.0024984486
## 
## Transition probability matrix (based on mean covariate values):
## ---------------------------------------------------------------
##           resident   transit
## resident 0.8839104 0.1160896
## transit  0.4393721 0.5606279
## 
## Initial distribution:
## ---------------------
##  resident   transit 
## 0.5030794 0.4969206

pr_high <- pseudoRes(hmm2_high_young_sens_bnsm)
qqnorm(pr_high$stepRes,
       pch = 21,          
       col = "midnightblue",   
       cex = 2.2,         
       main = "Q-Q Plot of Step Residuals",
       xlab = "Theoretical Quantiles",
       ylab = "Sample Quantiles",
       cex.main = 2.2,      
       cex.lab = 1.8,     
       cex.axis = 1.6)

qqnorm(pr_high$angleRes,
       pch = 21,          
       col = "midnightblue",   
       cex = 2.2,         
       main = "Q-Q Plot of Turning angles Residuals",
       xlab = "Theoretical Quantiles",
       ylab = "Sample Quantiles",
       cex.main = 2.2,      
       cex.lab = 1.8,     
       cex.axis = 1.6)

4.3.2 Stationary Plots (High)

plotStationary(
  hmm2_high_young_sens_bnsm, 
  plotCI = TRUE, 
  col = c("#0072B2", "#d33420"),  
  cex.axis = 2, 
  cex.lab = 2,
  cex.legend = 1.8,
  cex.main = 2,
  lwd =2,
  legend.pos = "topright")

4.4 Proportions of each state

low <- as.data.frame(low) %>% 
  mutate(state2 = as.factor(hmm_2states_low))
high <- as.data.frame(high) %>% 
  mutate(state2 = as.factor(hmm_2states_high))

color <- c("Native Forest Fragment" = "palegreen4", "Plantation Stand" = "saddlebrown")

state_proportions_low <- low %>%
  group_by(tipouso_grouped, state2) %>%
  summarise(count = n(), .groups = "drop") %>%
  group_by(tipouso_grouped) %>%
  mutate(proportion = count / sum(count)) %>%
  arrange(tipouso_grouped, desc(proportion))

state_proportions_low <- state_proportions_low %>%
  mutate(tipouso_grouped = recode(tipouso_grouped, 
                                  "BNSM" = "Native Forest Fragment",
                                  "plantations" = "Plantation Stand"))
state_proportions_low <- state_proportions_low %>%
  mutate(state2 = recode(state2, 
                                  "1" = "Resident",
                                  "2" = "Transit")) %>% 
  mutate(season = "Low Bark-stripping")

# High

state_proportions_high <- high %>%
  group_by(tipouso_grouped, state2) %>%
  summarise(count = n(), .groups = "drop") %>%
  group_by(tipouso_grouped) %>%
  mutate(proportion = count / sum(count)) %>%
  arrange(tipouso_grouped, desc(proportion))

state_proportions_high <- state_proportions_high %>%
  mutate(tipouso_grouped = recode(tipouso_grouped, 
                                  "BNSM" = "Native Forest Fragment",
                                  "plantations" = "Plantation Stand"))
state_proportions_high <- state_proportions_high %>%
  mutate(state2 = recode(state2, 
                                  "1" = "Resident",
                                  "2" = "Transit"))%>% 
  mutate(season = "High Bark-stripping")

divided_state_props <- rbind (state_proportions_low, state_proportions_high) %>% 
  mutate(season = factor(season, levels = c("Low Bark-stripping", "High Bark-stripping")))

ggplot(divided_state_props, aes(x = state2, y = proportion, fill = tipouso_grouped)) +
  geom_bar(stat = "identity", position = "dodge", color = "black", size = 1.1) + 
  scale_fill_manual(values = color, 
                    name = "Habitat Type",  
                    labels = c("Native Forest Fragment", "Plantation Stand")) + 
  facet_wrap(~season) +
  theme_bw() +
  theme(text = element_text(family = "Calibri", size = 30),
    panel.grid.major = element_blank(),  
    panel.grid.minor = element_blank(),  
    legend.position = "bottom",  
    legend.title = element_text(size = 24), 
    legend.text = element_text(size = 24),  
    legend.key.size = unit(1.1, "cm"),  
    axis.text.x = element_text(size = 22, color = "black"),
    axis.text.y = element_text(size = 22, color = "black"),
    axis.title.x = element_blank(),
    plot.margin = margin(1, 1, 1, 1, "cm"))

4.5 Step Lenght and Turning Angles

hmm_prep_nf <- hmm_prep %>% filter (tipouso_grouped == "BNSM")
hmm_prep_pita <- hmm_prep %>% filter (tipouso_grouped == "plantations")

sl_plant_nf <- ggplot() +
  
  geom_histogram(
    data = subset(hmm_prep_nf, state2 == "Resident"),
    aes(x = step, fill = state2),
    bins = 95, 
    alpha = 0.7,  
    color = "#013c5d",
    size=0.3,
    show.legend = FALSE) +
  
 
  geom_histogram(
    data = subset(hmm_prep_nf, state2 == "Transit"),
    aes(x = step, fill = state2),
    bins = 95, 
    alpha = 0.9,  
    color = "#6b170d",
    size=0.3,
    show.legend = FALSE) +
  
 
  scale_fill_manual(
    values = c(
      "Transit" = "#d33420", 
      "Resident" = "#0072B2")) +
  labs(
    x = "Step Length (m)",
    y = "Count") +
  scale_x_continuous(limits = c(0, 600)) +
  theme_minimal() +
  theme(
    text = element_text(family = "Calibri", size = 38, color = 'black'),
    panel.background = element_rect(fill = "white", color = NA),
    plot.background = element_rect(fill = "white", color = NA),
    axis.text = element_text(color = "black"),
    legend.text = element_text(size = 35),           
    legend.title = element_text(size = 35),         
    axis.title.x = element_text(vjust = 0.5, size =42),         
    axis.title.y = element_text(vjust = 2, size =42),       
    strip.text = element_text(size = 35),
    legend.position = "none", 
    plot.margin = unit(c(1.5, 1.5, 1.5, 1.5), "cm"))   

roset_nf <- ggplot() +
  geom_histogram(
    data = subset(hmm_prep_nf, state2 == "Resident"),
    aes(x = angle, fill = state2),
    bins = 45, 
    alpha = 0.7,  
    color = "#013c5d",
    size=0.3,
    show.legend = FALSE) +
  geom_histogram(
    data = subset(hmm_prep_nf, state2 == "Transit"),
    aes(x = angle, fill = state2),
    bins = 45, 
    alpha = 0.9,  
    color = "#6b170d",
    size=0.3,
    show.legend = FALSE) +

  scale_fill_manual(
    values = c(
      "Transit" = "#d33420",  
      "Resident" = "#0072B2")) +

  coord_polar(start = -pi) +
  

  labs(
    x = "Turning angle (radians)", 
    y = "Count") + 

  scale_x_continuous(
    breaks = c(-pi/2, 0, pi/2, pi),
    labels = c(
      expression(-frac(pi, 2)),  
      "0",
      expression(frac(pi, 2)),  
      expression(pi))) +

  theme_minimal() +
  theme(
    text = element_text(family = "Calibri", size = 42),
    panel.background = element_rect(fill = "white", color = NA),  
    plot.background = element_rect(fill = "white", color = NA),  
    panel.grid = element_line(color = "gray80"),  
    axis.text.x = element_text(size = 35, color = "black", hjust = 10), 
    axis.text.y = element_text(size = 35, color = "black"), 
    axis.title.x = element_text(size = 42, margin = margin(t = 35)), 
    axis.title.y = element_text(size = 42, margin = margin(r = 20)),
    legend.position = "none",
    plot.margin = unit(c(1.5, 1.5, 1.5, 1.5), "cm"))


sl_plant_pita <- ggplot() +
  
  geom_histogram(
    data = subset(hmm_prep_pita, state2 == "Resident"),
    aes(x = step, fill = state2),
    bins = 95, 
    alpha = 0.7,  
    color = "#013c5d",
    size=0.3,
    show.legend = FALSE) +
  
 
  geom_histogram(
    data = subset(hmm_prep_pita, state2 == "Transit"),
    aes(x = step, fill = state2),
    bins = 95, 
    alpha = 0.9,  
    color = "#6b170d",
    size=0.3,
    show.legend = FALSE) +
    scale_fill_manual(
    values = c(
      "Transit" = "#d33420", 
      "Resident" = "#0072B2")) +
  labs(
    x = "Step Length (m)",
    y = "Count") +
  scale_x_continuous(limits = c(0, 600)) +
  theme_minimal() +
   theme(
    text = element_text(family = "Calibri", size = 38, color = 'black'),
    panel.background = element_rect(fill = "white", color = NA),
    plot.background = element_rect(fill = "white", color = NA),
    axis.text = element_text(color = "black"),
    legend.text = element_text(size = 35),           
    legend.title = element_text(size = 35),         
    axis.title.x = element_text(vjust = 0.5, size =42),         
    axis.title.y = element_text(vjust = 2, size =42),       
    strip.text = element_text(size = 35),
    legend.position = "none", 
    plot.margin = unit(c(1.5, 1.5, 1.5, 1.5), "cm"))   

roset_pita <- ggplot() +
  geom_histogram(
    data = subset(hmm_prep_pita, state2 == "Transit"),
    aes(x = angle, fill = state2),
    bins = 45, 
    alpha = 0.9,  
    color = "#6b170d",
    size=0.3,
    show.legend = FALSE) +
 
  geom_histogram(
    data = subset(hmm_prep_pita, state2 == "Resident"),
    aes(x = angle, fill = state2),
    bins = 45, 
    alpha = 0.7,  
    color = "#013c5d",
    size=0.3,
    show.legend = FALSE) +

  scale_fill_manual(
    values = c(
      "Transit" = "#d33420",
      "Resident" = "#0072B2")) +

  coord_polar(start = -pi) +
  
  
  labs(
    x = "Turning angle (radians)", 
    y = "Count") + 

  scale_x_continuous(
    breaks = c(-pi/2, 0, pi/2, pi),
    labels = c(
      expression(-frac(pi, 2)),  
      "0",
      expression(frac(pi, 2)),  
      expression(pi))) +

  theme_minimal() +
  theme(
    text = element_text(family = "Calibri", size = 42),
    panel.background = element_rect(fill = "white", color = NA),  
    plot.background = element_rect(fill = "white", color = NA),  
    panel.grid = element_line(color = "gray80"),  
    axis.text.x = element_text(size = 35, color = "black", hjust = 10), 
    axis.text.y = element_text(size = 35, color = "black"), 
    axis.title.x = element_text(size = 42, margin = margin(t = 35)), 
    axis.title.y = element_text(size = 42, margin = margin(r = 20)),
    legend.position = "none",
    plot.margin = unit(c(1.5, 1.5, 1.5, 1.5), "cm"))


combined_plot_habitats <- sl_plant_nf + roset_nf + sl_plant_pita + roset_pita +
  plot_layout(guides = "collect")   

combined_plot_habitats <- combined_plot_habitats & 
  theme(
    legend.position = "bottom",
    legend.title = element_text(size = 34, face = "bold"),
    legend.text = element_text(size = 32),
    legend.key.size = unit(1.5, "cm"),  
    legend.spacing.y = unit(0.9, "cm"))

combined_plot_habitats

Figure 2. Distributions of step lengths and turning angles for two behavioral movement states of black capuchin monkeys in pine plantations of northeastern Argentina. The Resident state (blue) is characterized by short step lengths and more random turning angles, and was more common in native forest fragments (upper panels). The Transit state (red) is characterized by longer step lengths and more linear movements, and was more common in pine stands (lower panels).

5 Iguazu National Park Data

5.1 INP processed data

hmm_inp <- readRDS("files/hmm_pni_final.rds")

5.2 Preparing data for HMM

hmm_inp <- 
  hmm_inp %>% filter(!is.na(lon) & 
                        !is.na(lat)) 
hmm_sf_inp <-
  hmm_inp %>% 
  st_as_sf(
    coords = c(
      'lon',
      'lat'),
    crs = 4326) %>% 
  st_transform(crs=32721) 

hmm_sf_inp_new_coords <- 
  hmm_sf_inp %>% 
  st_coordinates()

X <- hmm_sf_inp_new_coords[,"X"]
Y <- hmm_sf_inp_new_coords[,"Y"]

hmm_utm_inp <- cbind(
  hmm_inp,
  X,
  Y)


hmm_prep_inp <- 
  prepData(
    hmm_utm_inp, type = "UTM",
    coordNames = c("X","Y")) 

hmm_prep_inp <- 
  hmm_prep_inp %>% filter(!is.na(step) & 
                        !is.na(angle))

5.3 Initial Parameters

stepMean0_2 <- c(m1 = 20, m2 = 300)
stepSD0_2 <- c(sd1 = 20, sd2 = 300)
zeromass_2 <- c(zm1 = 0.05, zm2 =0.05) 
angleCon0_2 <- c(rho1  = 0.1, rho2 = 0.5)

stateNames_2 <- c("resident","transit")

dist_2 <- list(step = "gamma", angle = "wrpcauchy")
Par0_2 <- list(step=c(stepMean0_2, stepSD0_2), angle = c(angleCon0_2))

5.4 Model Fitting

hmm_fit_inp <-
  fitHMM(
    data = hmm_prep_inp,
    nbStates = 2, dist = dist_2,
    Par0 = Par0_2, stateNames = stateNames_2,
    formula = ~ 1)

pr <- pseudoRes(hmm_fit_inp)
qqnorm(pr$stepRes,
       pch = 21,          
       col = "midnightblue",   
       cex = 2.2,         
       main = "Q-Q Plot of Step Residuals",
       xlab = "Theoretical Quantiles",
       ylab = "Sample Quantiles",
       cex.main = 2.2,      
       cex.lab = 1.8,     
       cex.axis = 1.6)

qqnorm(pr$angleRes,
       pch = 21,          
       col = "midnightblue",   
       cex = 2.2,         
       main = "Q-Q Plot of Turning angles Residuals",
       xlab = "Theoretical Quantiles",
       ylab = "Sample Quantiles",
     cex.main = 2.2,      
       cex.lab = 1.8,     
       cex.axis = 1.6)

hmm_fit_inp

## Value of the maximum log-likelihood: -19765.78 
## 
## 
## step parameters:
## ----------------
##      resident  transit
## mean 36.41865 97.31534
## sd   22.75771 59.20351
## 
## angle parameters:
## -----------------
##               resident   transit
## mean          0.000000 0.0000000
## concentration 0.238198 0.3410758
## 
## Regression coeffs for the transition probabilities:
## ---------------------------------------------------
##                1 -> 2    2 -> 1
## (Intercept) -3.388103 -1.813249
## 
## Transition probability matrix:
## ------------------------------
##           resident    transit
## resident 0.9673306 0.03266935
## transit  0.1402459 0.85975415
## 
## Initial distribution:
## ---------------------
##  resident   transit 
## 0.6928904 0.3071096

5.4.1 Proportions of each state

hmm_states_2_inp <- viterbi(hmm_fit_inp)


hmm_states_df_inp <- hmm_prep_inp %>% 
  mutate(state = hmm_states_2_inp) 


prop.table(table(hmm_states_2_inp))

## hmm_states_2_inp
##          1          2 
## 0.92508834 0.07491166

hmm_states_df_inp$state <- factor(
  hmm_states_df_inp$state,
  levels = c(2, 1),  # Resident first, Transit second
  labels = c("Transit", "Resident"))


sl_pni <- ggplot() +
  
  geom_histogram(
    data = subset(hmm_states_df_inp, state == "Resident"),
    aes(x = step, fill = state),
    bins = 95, 
    alpha = 0.7,  
    color = "#013c5d",
    size=0.3,
    show.legend = FALSE) +
  

  geom_histogram(
    data = subset(hmm_states_df_inp, state == "Transit"),
    aes(x = step, fill = state),
    bins = 95, 
    alpha = 0.9,  
    color = "#6b170d",
    size=0.3,
    show.legend = FALSE) +
  

  scale_fill_manual(
    values = c(
      "Transit" = "#d33420",  
      "Resident" = "#0072B2")) +
  labs(
    x = "Step Length (m)",
    y = "Count") +
  scale_x_continuous(limits = c(0, 400)) +
  theme_minimal() +
 theme(
    text = element_text(family = "Calibri", size = 38, color = 'black'),
    panel.background = element_rect(fill = "white", color = NA),
    plot.background = element_rect(fill = "white", color = NA),
    axis.text = element_text(color = "black"),
    legend.text = element_text(size = 35),           
    legend.title = element_text(size = 35),         
    axis.title.x = element_text(vjust = 0.5, size =42),         
    axis.title.y = element_text(vjust = 2, size =42),       
    strip.text = element_text(size = 35),
    legend.position = "none", 
    plot.margin = unit(c(1.5, 1.5, 1.5, 1.5), "cm")) 

roset_pni <- ggplot() +
  geom_histogram(
    data = subset(hmm_states_df_inp, state == "Resident"),
    aes(x = angle, fill = state),
    bins = 45, 
    alpha = 0.7,  
    color = "#013c5d",
    size=0.3,
    show.legend = FALSE) +
  geom_histogram(
    data = subset(hmm_states_df_inp, state == "Transit"),
    aes(x = angle, fill = state),
    bins = 45, 
    alpha = 0.9,  
    color = "#6b170d",
    size=0.3,
    show.legend = FALSE) +

  scale_fill_manual(
    values = c(
      "Transit" = "#d33420",  
      "Resident" = "#0072B2")) +

  coord_polar(start = -pi) +
  
  # Labels
  labs(
    x = "Turning angle (radians)", 
    y = "Count") + 

  scale_x_continuous(
    breaks = c(-pi/2, 0, pi/2, pi),
    labels = c(
      expression(-frac(pi, 2)),  
      "0",
      expression(frac(pi, 2)),  
      expression(pi))) +

  theme_minimal() +
  theme(
    text = element_text(family = "Calibri", size = 42),
    panel.background = element_rect(fill = "white", color = NA),  
    plot.background = element_rect(fill = "white", color = NA),  
    panel.grid = element_line(color = "gray80"),  
    axis.text.x = element_text(size = 35, color = "black", hjust = 10), 
    axis.text.y = element_text(size = 35, color = "black"), 
    axis.title.x = element_text(size = 42, margin = margin(t = 35)), 
    axis.title.y = element_text(size = 42, margin = margin(r = 20)),
    legend.position = "none",
    plot.margin = unit(c(1.5, 1.5, 1.5, 1.5), "cm"))

combined_plot_inp <- sl_pni + roset_pni +
  plot_layout(guides = "collect")   

combined_plot_inp <- combined_plot_inp +
  theme(
    plot.title = element_text(
      size = 39,
      hjust = 1,
      margin = margin(b = 50)),
    legend.position = "bottom",
    legend.title = element_text(size = 34, face = "bold"),
    legend.text = element_text(size = 32),
    legend.key.size = unit(1.5, "cm"),
    legend.spacing.y = unit(0.9, "cm"))
combined_plot_inp

HMM Iguazu National Park. Distributions of step lengths and turning angles for two behavioral movement states of black capuchin monkeys in Iguazu National Park, northeastern Argentina. The Resident state (blue) is characterized by short step lengths and more random turning angles. The Transit state (red) is characterized by longer step lengths and more linear movements.