We compared the daily trajectories of black capuchins in pine plantations vs INP to investigate how travel distances reflect movement patterns and ecological requirements in contrasting landscapes. The availability and spatial distribution of sleeping sites are recognized as critical factors influencing the daily movement patterns of arboreal primates, determining where individuals begin and end their activity periods (Fruth et al., 2018). To integrate this factor into our analysis, we calculated the straight-line distance between the morning and evening sleeping sites, which were defined as locations recorded within 15 min after sunrise or before sunset. Daily sunrise and sunset times were computed using the ‘suncalc’ R package (Thieurmel et al., 2019), ensuring temporal accuracy across the study period. Only complete tracking days were included in the analysis, defined as days with eight evenly spaced locations recorded at 90-minute intervals (±15 minutes), capturing the full extent of daily activity.

1 R Packages

library(amt)
library(tidyverse)
library(sf)
library(adehabitatLT)
library(fitdistrplus)
library(cowplot)
library(lmtest)
library(MuMIn)
library(lme4)
library(caret)
library(flextable)
library(glmmTMB)
library(purrr)
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")
showtext_auto()

2 Raw DataSet

locs <- readRDS('files/hr.rds')%>% 
  rename (id = individual.local.identifier,
          lat = location.lat,
          lon = location.long,
          t = timestamp)

locs <- locs %>%
  mutate(
    t = as.POSIXct(sprintf("%s-%02d-%02d %02d:%02d:%02d", 
                           year, 
                           month(t), 
                           day(t), 
                           hour(t), 
                           minute(t), 
                           second(t)),
                   tz = attr(t, "tzone")))

locs %>% group_by(area) %>% summarise(locations = n())

Including the Shape file of Timber companies to have landscape features

shp <- st_read("files/UNION.shp")

crs_shp <- st_crs(shp)

locs_ara <- locs %>% filter (area == "Arauco")
locs_pni <- locs %>% filter (area == "PN Iguazu")


locs_ara_sf <- locs_ara %>%  
  st_as_sf(coords = c('lon', 'lat'),
           crs = 4326) %>% 
  st_transform(crs_shp)

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


locs_ara_sf <- locs_ara_sf %>%
  st_join(shp[, c("tipouso_grouped")])

locs_ara_sf <- st_transform(locs_ara_sf, 4326)

locs_ara2 <- locs_ara_sf %>%
  mutate(lon = st_coordinates(.)[, 1],
         lat = st_coordinates(.)[, 2]) %>%
  st_drop_geometry()

Unifying the processed data, and make Track XYT object

locs <- bind_rows(locs_ara2,locs_pni)

locs_track <-  locs %>%
  make_track(.x = lat,
             .y = lon,
             .t = t,
             id = id,
             year = year,
             grupo = grupo,
             month = month,
             tipouso_grouped = tipouso_grouped,
             area=area,
             crs = 32721,
             all_cols = TRUE)

Continue with data processing

locs <- locs %>%
  mutate(date = as.Date(t)) %>%  
  group_by(id, year, date) %>%
  mutate(count_locs = n()) %>%  
  ungroup()


locs <- locs %>%
  mutate(id = case_when(
    id == "spo" & year == 2022 ~ "spo.22",
    id == "spo" & year == 2014 ~ "spo.14",  
    TRUE ~ id))


locs <- locs %>%
  mutate(
    t = case_when(
      id %in% c("col.oto", "gua.tilo", "pan.pri", "tri.pri", 'spo.22') ~ t - lubridate::hours(3),
      TRUE ~ t ))

locs_clean <- locs %>%
  filter(count_locs > 8)


intervals <- locs_clean %>%
  group_by(id, year, date) %>%  
  arrange(t, .by_group = TRUE) %>% 
  mutate(
    interval_minutes = as.numeric(difftime(t, lag(t), units = "mins"))
  ) %>%
  ungroup()




intervals_filtered <- intervals %>% filter(!is.na(interval_minutes))


locs_clean_filtered <- intervals %>%
  filter(is.na(interval_minutes) | interval_minutes <= 200) %>%  
  group_by(id, year, date) %>%
  filter(lead(interval_minutes, default = 0) <= 200 | is.na(lead(interval_minutes))) %>%  
  ungroup()

locs_clean_filtered <- locs_clean_filtered %>%
  mutate(
    id2 = id,  
    id = paste0(tolower(substr(grupo, 1, 3)), ".",  
                tolower(substr(month, 1, 3)), ".", 
                substr(year, 3, 4)))

locs_clean_filtered <- 
  locs_clean_filtered %>%
  mutate(date = as.Date(t)) %>% 
  group_by(id, date) %>%
  arrange(t, .by_group = TRUE) %>%  
  mutate(
    daily_track_min = as.numeric(difftime(max(t), min(t), units = "mins")), 
    daily_track_h = as.numeric(difftime(max(t), min(t), units = "hours"))    
  ) %>%
  ungroup()


locs_full_day <- locs_clean_filtered %>%
  filter(daily_track_h > 8)

date_counts <- locs_full_day %>%
  group_by(area, id) %>%
  summarise(num_dates = n_distinct(date))

total_dates_by_area <- date_counts %>%
  group_by(area) %>%
  summarise(total_num_dates = sum(num_dates), .groups = "drop")

locs_clean_filtered <- locs_clean_filtered %>%
  mutate(area = recode(area, "Arauco" = "Plantations", "PN Iguazu" = "INP"))
ggplot(intervals_filtered, aes(x = interval_minutes)) +
  geom_histogram(binwidth = 8, fill = "midnightblue", color = "black", alpha = 0.7) +  
  labs(
    title = "Intervals between locations distribution",
    x = "Interval (min)",
    y = "Frequency") +
  theme_minimal() +
  theme(text = element_text(family = 'Calibri', size =34))

we still have some big gaps in the data, intervals up to 400 min.

locs_full_day <- locs_clean_filtered %>%
  filter(daily_track_h > 8)

date_counts <- locs_full_day %>%
  group_by(area, id) %>%
  summarise(num_dates = n_distinct(date))

total_dates_by_area <- date_counts %>%
  group_by(area) %>%
  summarise(total_num_dates = sum(num_dates), .groups = "drop")

2.1 Calculating Steps

First I am creating the SF object to extract UTM coordinates (\(X\) and \(Y\)), using complex number then I obtain \(Z = X + iY\). The total daily trajectory includes the sum of all steps conducted within a same day. There are still some strange data that showed length above 7000 m, these days were removed becuase there were linked to location that not are from the monkeys (inside buildings for instance)

locs_sf <- locs_full_day %>% 
  st_as_sf(coords = c('lon', 'lat'), crs = 4326) %>% 
  st_transform(crs = 32721)

locs_sf_coords <- st_coordinates(locs_sf)

locs_full_day <- locs_full_day %>% 
  mutate(
    X = locs_sf_coords[,"X"],
    Y = locs_sf_coords[,"Y"],
    Z = X + 1i * Y
  )

steps_data <- locs_full_day %>%
  group_by(area, id, date) %>%
  arrange(t, .by_group = TRUE) %>%  
  mutate(S = c(Mod(diff(Z)), NA)) %>%  
  summarise(
    mean_S = mean(S, na.rm = TRUE),   
    sd_S = sd(S, na.rm = TRUE),       
    sum_S = sum(S, na.rm = TRUE),
    prop_plantations = mean(tipouso_grouped == "plantations", na.rm = TRUE),
    .groups = "drop"
  )

steps_data <- steps_data %>%
  filter(sum_S < 7000)

Exploring, little processing is still needed

steps_data <- steps_data %>%
  mutate(area = recode(area, "Arauco" = "Plantations", "PN Iguazu" = "INP"))

ggplot(steps_data, aes(x = area, y = sum_S, fill = area)) +
  geom_jitter(alpha = 0.1, width = 0.2, size= 3, height = 20) +
  geom_boxplot(alpha = 0.7, lwd = 1.2, outlier.shape = NA) +  
  scale_fill_manual(values = c("Plantations" = "#8B4513", "INP" = "#228B22")) +
  scale_y_continuous(limits = c(0, 3700), breaks = seq(0, 3700, by = 300)) +
  labs(
    x = "Study Site",
    y = "Daily trajectory (m)") +
  theme_minimal() +
  theme(legend.position = "none",
        text = element_text(family = 'Calibri', size = 34),
        axis.title.x = element_text(vjust = -4),
        axis.title.y = element_text(vjust = 5),
        panel.background = element_rect(fill = "white", color = NA),
        plot.background = element_rect(fill = "white", color = NA),
        plot.margin = unit(c(1, 1, 1, 1), "cm"))

2.2 Speeds

locs_full_day <- locs_full_day %>%
  mutate(
    speed = ifelse(!is.na(S) & !is.na(interval_minutes) & interval_minutes > 0, S / interval_minutes, NA))
speed <- locs_full_day %>%
  group_by(area, id, date) %>%
  arrange(t, .by_group = TRUE) %>%
  summarise(
    mean_S = mean(speed, na.rm = TRUE),   
    sd_S = sd(S, na.rm = TRUE), 
    prop_plantations = mean(tipouso_grouped == "plantations", na.rm = TRUE),
    .groups = "drop")
speed %>% 
  filter(area=='Plantations') %>% 
  ggplot(aes(x = prop_plantations, y = speed)) +
  geom_jitter(width= 0.01, alpha = 0.5, size = 2, color = "grey24") +  # Puntos para cada día
  geom_smooth(method = "glm", color = "darkred", size = 3, se = TRUE) +
  labs(x = "Daily proportions of locacions within plantations",
    y = "Mean daily speed (m/min)"
  ) +
  theme_minimal() +
  theme(legend.position = "none",
        text = element_text(family = 'Calibri', size = 34),
        axis.title.x = element_text(vjust = -4),
        axis.title.y = element_text(vjust = 5),
        panel.background = element_rect(fill = "white", color = NA),
        plot.background = element_rect(fill = "white", color = NA),
        plot.margin = unit(c(1, 1, 1, 1), "cm"))

2.3 Suncalc estimations

locs_first_last <- locs_full_day %>%
  group_by(id, date) %>%
  slice_min(order_by = t, n = 1) %>%
  mutate(time_label = "first") %>%  
  bind_rows(
    locs_full_day %>%
      group_by(id, date) %>%
      slice_max(order_by = t, n = 1) %>%
      mutate(time_label = "last")   
  ) %>%
  ungroup()

locs_full_day <- locs_full_day %>%
  mutate(date = as.Date(date))

get_sun_times <- function(date, lat, lon) {
  sun_times <- getSunlightTimes(
    date = date, lat = lat, lon = lon,
    keep = c("sunrise", "sunset"),
    tz = "America/Argentina/Buenos_Aires")
  return(sun_times)
}

sun_times_list <- pmap(
  list(locs_full_day$date, locs_full_day$lat, locs_full_day$lon),
  get_sun_times)

sun_data <- bind_rows(sun_times_list)

locs_full_day <- bind_cols(locs_full_day, sun_times_df)

locs_full_day <- locs_full_day %>%
  mutate(
    sunrise = format(sunrise, "%H:%M:%S"),
    sunset = format(sunset, "%H:%M:%S"))

locs_first_last <- locs_first_last %>%
  mutate(month_join = month(t))

sun_data <- sun_data %>%
  mutate(month_join = month(date))

locs_first_last <- locs_first_last %>%
  mutate(
    horario = format(t, format = "%H:%M:%S"))

sun_data <- sun_data %>%
  mutate(
    sunrise_horario = format(sunrise_datetime, format = "%H:%M:%S"),
    sunset_horario = format(sunset_datetime, format = "%H:%M:%S"))

sun_data <- sun_data %>%
  mutate(month_join = month(date)) %>%
  group_by(month_join) %>%
  slice(1) %>%  
  ungroup() %>%
  mutate(
    sunrise_horario = format(sunrise_datetime, format = "%H:%M:%S"),
    sunset_horario = format(sunset_datetime, format = "%H:%M:%S"))


locs_with_sun_info <- 
  locs_first_last %>%
  left_join(sun_data_unique, by = "month_join") %>%
  mutate(
    
    t_horario = as.POSIXct(paste("2024-01-01", horario)),
    sunrise_period = as.POSIXct(paste("2024-01-01", sunrise_horario)),
    sunset_period = as.POSIXct(paste("2024-01-01", sunset_horario)),
    
    t_to_sunrise = if_else(
      time_label == "first",
      as.numeric(difftime(
        t_horario,
        sunrise_period,
        units = "mins")),
      NA_real_),
    t_to_sunset = if_else(
      time_label == "last", 
      as.numeric(difftime(
        sunset_period,
        t_horario, 
        units = "mins")),
      NA_real_)) 

locs_with_sun_info <- locs_with_sun_info %>%
  dplyr::select(-sunrise_period, sunset_period, sunset,sunrise)

locs_with_sun_info <- locs_with_sun_info %>%
  mutate(
    time_to_sun = coalesce(t_to_sunrise, t_to_sunset)) %>%
  dplyr::select(-t_to_sunrise, -t_to_sunset)


locs_with_sun <- locs_with_sun_info %>%
  filter(time_to_sun >= -15 & time_to_sun <= 15)

2.4 Distances between sleeping sites (euclicean distances)

locs_with_sun <- locs_with_sun %>%
  group_by(id, date.x) %>%
  mutate(
    X_first = na_if(first(X[time_label == "first"]), NA),
    Y_first = na_if(first(Y[time_label == "first"]), NA),
    X_last = na_if(first(X[time_label == "last"]), NA),
    Y_last = na_if(first(Y[time_label == "last"]), NA),
    
    distance_ss = if_else(
      !is.na(X_first) & !is.na(Y_first) & !is.na(X_last) & !is.na(Y_last),
      sqrt((X_last - X_first)^2 + (Y_last - Y_first)^2),
      NA_real_
    )
  ) %>%
  ungroup() %>%
  dplyr::select(-X_first, -Y_first, -X_last, -Y_last)

Checking some specific data

demos1 <- locs_full_day %>%
  filter((id == "gua.may.22" & date == as.Date("2022-05-25")))%>%
  mutate(hour = hour(t)) 
           
demos2 <- locs_full_day %>%
           filter((id == "gun.nov.10" & date == as.Date("2010-11-22")))%>%
  mutate(hour = hour(t)) 


sf_demos1 <- st_as_sf(demos1, coords = c("lon", "lat"), crs = 4326)
sf_demos2 <- st_as_sf(demos2, coords = c("lon", "lat"), crs = 4326)


first_point1 <- sf_demos1 %>% filter(t == min(t))
last_point1 <- sf_demos1 %>% filter(t == max(t))

first_point1_utm <- st_transform(first_point1, crs = 32721)  
last_point1_utm <- st_transform(last_point1, crs = 32721)

distance_meters1 <- st_distance(first_point1_utm, last_point1_utm)
distance_meters1 <- as.numeric(distance_meters1)

mid_x <- (st_coordinates(first_point1)[1,1] + st_coordinates(last_point1)[1,1]) / 2
mid_y <- (st_coordinates(first_point1)[1,2] + st_coordinates(last_point1)[1,2]) / 2



first_point2 <- sf_demos2 %>% filter(t == min(t))
last_point2 <- sf_demos2 %>% filter(t == max(t))

first_point2_utm <- st_transform(first_point2, crs = 32721)  
last_point2_utm <- st_transform(last_point2, crs = 32721)

distance_meters2 <- st_distance(first_point2_utm, last_point2_utm)
distance_meters2 <- as.numeric(distance_meters2)

mid_x2 <- (st_coordinates(first_point2)[1,1] + st_coordinates(last_point2)[1,1]) / 2
mid_y2 <- (st_coordinates(first_point2)[1,2] + st_coordinates(last_point2)[1,2]) / 2

sf_demos1_df <- sf_demos1 %>%
  mutate(lon = st_coordinates(.)[,1], lat = st_coordinates(.)[,2])

sf_demos2_df <- sf_demos2 %>%
  mutate(lon = st_coordinates(.)[,1], lat = st_coordinates(.)[,2])
ggplot() +
  geom_path(data = sf_demos1_df, aes(x = lon, y = lat), color = "grey23", size = 1.4, linetype = "solid") +
  geom_sf(data = sf_demos1, size = 4, color = "grey5") +  
  geom_segment(data = NULL,
               aes(x = st_coordinates(first_point1)[1,1],
                   y = st_coordinates(first_point1)[1,2],
                   xend = st_coordinates(last_point1)[1,1],
                   yend = st_coordinates(last_point1)[1,2]),
               color = "midnightblue", linetype = "dashed", size = 1.4) +
  annotate("text", x = mid_x, y = mid_y + 0.0002, label = paste0(round(distance_meters1, 1), " m"),
           color = "midnightblue", size = 8, angle = -12.5, fontface = "bold") +
  annotate("text", x = st_coordinates(first_point1)[1,1] + 0.0008, 
           y = st_coordinates(first_point1)[1,2] + 0.00005, 
           label = "Sunrise", color = "goldenrod4", size = 8, fontface = "bold") +
  annotate("text", x = st_coordinates(last_point1)[1,1], 
           y = st_coordinates(last_point1)[1,2] - 0.00025, 
           label = "Sunset", color = "darkred", size = 8, fontface = "bold") +
  labs(title = "Daily trakectory of Columba, 2022-04-24",
       x = "Longitude",
       y = "Latitud") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, vjust = 0.5),
        legend.position = "top",
        legend.title = element_blank(),
        text = element_text(family = 'Calibri', size = 34),
        plot.title = element_text(hjust = 0.5),
        axis.title.x = element_text(vjust = -3),
        axis.title.y = element_text(vjust = 4),
        panel.background = element_rect(fill = "white", color = NA),
        plot.background = element_rect(fill = "white", color = NA),
        plot.margin = unit(c(1, 1, 1, 1), "cm"))

We will keep only complete tracking days, defined as days with eight evenly spaced locations recorded at 90-minute intervals (±15 minutes).

locs_day <- locs_with_sun %>%
  group_by(id, date) %>%
  mutate(
    time_diff = difftime(t, lag(t), units = "mins"),
    valid_intervals = sum(abs(time_diff - 90) <= 15, na.rm = TRUE)) %>%
  filter(n() == 8 & valid_intervals == 7) %>%
  ungroup()
ggplot(locs_with_sun, aes(x = distance_ss, fill = area)) +
  geom_density(alpha = 0.5) +
  scale_fill_manual(values = c("Plantations" = "#8B4513", "INP" = "#228B22")) + 
  scale_y_continuous(labels = scales::comma)+
  scale_x_continuous(limits = c(0,3000))+
  labs(
    x = "Distance from morning to evening sleeping site (m)",
    y = "Density")+
    guides(fill = guide_legend(title = NULL)) +
  theme_minimal()+
  theme(text = element_text(family = 'Calibri', size = 34),
        axis.title.x = element_text(vjust = -4),
        axis.title.y = element_text(vjust = 5),
        panel.background = element_rect(fill = "white", color = NA),
        plot.background = element_rect(fill = "white", color = NA),
        plot.margin = unit(c(1, 1, 1, 1), "cm"))

2.5 Wheather Covariates

Minimum, mean and maximum daily temperature and precipitation (obtained for the Iguazú International Airport by the National Meteorological Service of Argentina), as these climatic factors may influence primate movement patterns (Asensio et al., 2012)

snm <- read_csv('igua_aero.csv') 

snm$date <- as.Date(snm$date, format = "%d/%m/%Y")

locs_day <- locs_day %>%
  left_join(snm, by = "date")

3 GLMM models

colores <- c('Plantations' = 'orange3', 'INP' = 'palegreen4')
colores2 <- c('Plantations' = '#f3c554', 'INP' = '#afeeb3')

3.1 Response variable distribution

normal_fit <- fitdist(locs_day$sum_S2, "norm")
lognormal_fit <- fitdist(locs_day$sum_S2, "lnorm")
gamma_fit <- fitdist(locs_day$sum_S2, "gamma")

cdf_sum_S2 <- cdfcomp(
  list(lognormal_fit, normal_fit, gamma_fit),
  main = "",
  legendtext = c("Log-normal", "Normal", "Gamma"),
  fitcol = c("darkred", "orange3", "goldenrod1"),
  plotstyle = "ggplot") +
  geom_line(size = 1.2) +
  theme(
    axis.title = element_text(size = 42),
    axis.text = element_text(size = 42),
    axis.text.x = element_text(angle = 45, vjust = 0.8),
    legend.position = c(0.8, 0.8),
    legend.text = element_text(size = 40))

qq_sum_S2 <- qqcomp(
  list(lognormal_fit, normal_fit, gamma_fit),
  main = "",
  fitcol = c("darkred", "orange3", "goldenrod1"),
  plotstyle = "ggplot") +
  geom_line(size = 1.2) +
  theme_bw() +
  theme(
    axis.title = element_text(size = 42),
    axis.text = element_text(size = 42),
    axis.text.x = element_text(angle = 45, vjust = 0.8),
    legend.position = "none")

plot_grid(qq_sum_S2, cdf_sum_S2, NULL, ncol = 3, rel_widths = c(1.5, 1.5, 0.75))

aic_values <- data.frame(
  Distribution = c("Normal", "Lognormal", "Gamma"),
  AIC = c(normal_fit$aic, lognormal_fit$aic, gamma_fit$aic)) %>%
  mutate(deltaAIC = AIC - min(AIC)) %>%
  arrange(AIC)

aic_table <- aic_values %>%
  flextable() %>%
  theme_tron_legacy() %>%  
  colformat_num(j = c("AIC", "deltaAIC"), digits = 3) %>%
  bold(part = "header") %>%
  autofit() 

aic_table

Distribution

AIC

deltaAIC

Normal

9,726.277

0.00000

Gamma

9,786.895

60.61803

Lognormal

9,848.748

122.47118

AIC table Best fit distributions for response variable (daily travel)

locs_day <- locs_day %>%
  mutate(
    distance_ss_scaled = scale(distance_ss, center = TRUE, scale = TRUE),
    t_mean_scaled = scale(t_mean, center = TRUE, scale = TRUE),
    precip_scaled = scale(precip, center = TRUE, scale = TRUE),
    t_min_scaled = scale(t_min, center = TRUE, scale = TRUE))

locs_day <- locs_day %>% 
  mutate(food = case_when(
    estacion2 == "invierno" ~ "low",
    TRUE ~ "high"))

3.2 Models

m0 <- lmer(sum_S2 ~ 1 + (1 | grupo), data = locs_day)
m1 <- lmer(sum_S2 ~ area + (1 | grupo), data = locs_day)
m2 <- lmer(sum_S2 ~ area + distance_ss_scaled + (1 | grupo), data = locs_day)
m3 <- lmer(sum_S2 ~ area * distance_ss_scaled + (1 | grupo), data = locs_day)
m4 <- lmer(sum_S2 ~ area * food + distance_ss_scaled +  (1 | grupo), data = locs_day)
m5 <- lmer(sum_S2 ~ area * distance_ss_scaled + food +  (1 | grupo), data = locs_day)
m6 <- lmer(sum_S2 ~ area + t_min_scaled + (1 | grupo), data = locs_day)
m7 <- lmer(sum_S2 ~ area + t_min_scaled + precip_scaled + (1 | grupo), data = locs_day)
m8 <- lmer(sum_S2 ~ area + distance_ss_scaled + t_min_scaled + precip_scaled + (1 | grupo), data = locs_day)
m9 <- lmer(sum_S2 ~ area * distance_ss_scaled + t_min_scaled + precip_scaled + (1 | grupo), data = locs_day)
m10 <- lmer(sum_S2 ~ area  + precip_scaled + distance_ss_scaled + t_min_scaled + (1 | grupo), data = locs_day)
m11 <- lmer(sum_S2 ~ area + t_min_scaled * distance_ss_scaled + precip_scaled + (1 | grupo), data = locs_day)
m12 <- lmer(sum_S2 ~ area * distance_ss_scaled + food + t_min_scaled + precip_scaled + (1 | grupo), data = locs_day)
m13 <- lmer(sum_S2 ~ area + distance_ss_scaled + food + t_min_scaled + precip_scaled + (1 | grupo), data = locs_day)
m14 <- lmer(sum_S2 ~ area * distance_ss_scaled + precip_scaled + (1 | grupo), data = locs_day)
m15 <- lmer(sum_S2 ~ area * food + t_min_scaled + (1 | grupo), data = locs_day)
m16 <- lmer(sum_S2 ~ area + distance_ss_scaled + food + daily_length + precip_scaled + (1 | grupo), data = locs_day)
m17 <- lmer(sum_S2 ~ area * distance_ss_scaled + daily_length + (1 | grupo), data = locs_day)
m18 <- lmer(sum_S2 ~ area * food + daily_length + t_min_scaled + (1 | grupo), data = locs_day)
m19 <- lmer(sum_S2 ~ area * daily_length + distance_ss_scaled + food + precip_scaled + (1 | grupo), data = locs_day)
m20 <- lmer(sum_S2 ~ area * daily_length * distance_ss_scaled + food + precip_scaled + (1 | grupo), data = locs_day)
m21<- lmer(sum_S2 ~ area * distance_ss_scaled + daily_length + food + t_min_scaled + precip_scaled + (1 | grupo), data = locs_day)

3.3 AIC comparisons

model_aic <- AIC(m0, m1, m2, m3, m4, m5, m6, m7, m8, m9, m10, m11, m12, m13, m14, m15, m16, m17, m18, m19, m20, m21)

model_aic$deltaAIC <- model_aic$AIC - min(model_aic$AIC)

model_aic$Rank <- rank(model_aic$deltaAIC, ties.method = "first")

model_aic$Model <- c(
  "Null model", #m0
  "Study site", #m1
  "Study site + Sleeping Sites Distance", #m2
  "Study site * Sleeping Sites Distance", #m3
  "Study site + Fruit Season + Sleeping Sites Distance", #m4
  "Study site * Sleeping Sites Distance + Fruit Season", #m5
  "Study site + Min. Temperature", #m6
  "Study site + Min. Temperature + Precipitation", #m7
  "Study site + Sleeping Sites Distance + Min. Temperature + Precipitation", #m8
  "Study site * Sleeping Sites Distance + Min. Temperature + Precipitation", #m9
  "Study site + Sleeping Sites Distance + Min. Temperature + Precipitation", #m10
  "Study site + Min. Temperature * Sleeping Sites Distance + Precipitation", #m11 
  "Study site * Sleeping Sites Distance + Fruit Season + Min. Temperature + Precipitation", #m12
  "Study site + Sleeping Sites Distance + Fruit Season + Min. Temperature + Precipitation", #m13
  "Study site * Sleeping Sites Distance + Precipitation", #m14
  "Study site * Fruit Season + Min. Temperature", #m15
  "Study site + Sleeping Sites Distance + Fruit Season + Daily Length + Precipitation", #m16
  "Study site * Sleeping Sites Distance + Daily Length", #m17
  "Study site * Fruit Season + Daily lenght + Min. Temperature", #m18
  "Study site * Daily length + Sleeping Sites Distance + Fruit Season + Precipitation", #m19
  "Study site * Daily length * Sleeping Sites Distance + Fruit Season + Precipitation", #m20
  "Study site * Sleeping Sites Distance + Daily Length + Fruit Season + Min. Temperature + Precipitation") #m21

model_table <- model_aic[order(model_aic$Rank), c("Rank", "Model", "df", "AIC", "deltaAIC")]
model_table <- model_table %>% arrange(Rank)

model_table %>%
  flextable() %>%
  theme_tron_legacy() %>%
  colformat_num(j = c("AIC", "deltaAIC"), digits = 3) %>%
  bold(part = "header") %>%
  set_caption("AIC Model Comparison for Study Site Effects") %>%
  autofit()
AIC Model Comparison for Study Site Effects

Rank

Model

df

AIC

deltaAIC

1

Study site * Sleeping Sites Distance + Daily Length + Fruit Season + Min. Temperature + Precipitation

10

9,467.830

0.000000

2

Study site * Sleeping Sites Distance + Fruit Season + Min. Temperature + Precipitation

9

9,475.603

7.772724

3

Study site * Sleeping Sites Distance + Min. Temperature + Precipitation

8

9,488.406

20.576256

4

Study site * Fruit Season + Daily lenght + Min. Temperature

8

9,503.292

35.462014

5

Study site + Sleeping Sites Distance + Fruit Season + Min. Temperature + Precipitation

8

9,508.242

40.412077

6

Study site * Fruit Season + Min. Temperature

7

9,511.325

43.494519

7

Study site + Min. Temperature * Sleeping Sites Distance + Precipitation

8

9,513.471

45.641324

8

Study site + Sleeping Sites Distance + Min. Temperature + Precipitation

7

9,520.213

52.383413

9

Study site + Sleeping Sites Distance + Min. Temperature + Precipitation

7

9,520.213

52.383413

10

Study site + Min. Temperature + Precipitation

6

9,527.070

59.239589

11

Study site + Min. Temperature

5

9,533.672

65.842111

12

Study site * Daily length * Sleeping Sites Distance + Fruit Season + Precipitation

12

9,550.382

82.552308

13

Study site * Sleeping Sites Distance + Fruit Season

7

9,595.433

127.603250

14

Study site * Sleeping Sites Distance + Daily Length

7

9,599.418

131.588258

15

Study site * Daily length + Sleeping Sites Distance + Fruit Season + Precipitation

9

9,600.521

132.691392

16

Study site * Sleeping Sites Distance + Precipitation

7

9,600.847

133.016548

17

Study site * Sleeping Sites Distance

6

9,608.300

140.469961

18

Study site + Sleeping Sites Distance + Fruit Season + Daily Length + Precipitation

8

9,613.366

145.535870

19

Study site + Fruit Season + Sleeping Sites Distance

7

9,617.597

149.766521

20

Study site + Sleeping Sites Distance

5

9,640.340

172.510322

21

Study site

4

9,647.128

179.298477

22

Null model

3

9,666.865

199.035172

Table A.1. The table presents a set of candidate models ranked according to Akaike’s Information Criterion (AIC), with the best-supported model (lowest AIC) listed first. The response variable is the daily trajectory length of black capuchin monkeys (Sapajus nigritus), modeled as a function of different ecological and environmental covariates. The best-ranked model includes interactions between Study Site and Sleeping Sites, as well as Fruit Season, Minimum Temperature, and Precipitation. The ΔAIC values indicate the relative support for each model, and “df” the degrees of freedom.

3.4 Best Model summary

summary(m21)
## Linear mixed model fit by REML ['lmerMod']
## Formula: 
## sum_S2 ~ area * distance_ss_scaled + daily_length + food + t_min_scaled +  
##     precip_scaled + (1 | grupo)
##    Data: locs_day
## 
## REML criterion at convergence: 9447.8
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.88224 -0.70006  0.09481  0.64915  2.67262 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  grupo    (Intercept)   3839    61.96  
##  Residual             473853   688.37  
## Number of obs: 599, groups:  grupo, 8
## 
## Fixed effects:
##                                    Estimate Std. Error t value
## (Intercept)                        2283.761    635.102   3.596
## areaPlantations                     403.969     86.748   4.657
## distance_ss_scaled                  140.085     40.217   3.483
## daily_length                        -17.739     49.644  -0.357
## foodlow                             158.756     89.674   1.770
## t_min_scaled                          6.509     35.733   0.182
## precip_scaled                        35.505     29.682   1.196
## areaPlantations:distance_ss_scaled -318.562     63.795  -4.993
## 
## Correlation of Fixed Effects:
##             (Intr) arPlnt dstn__ dly_ln foodlw t_mn_s prcp_s
## areaPlnttns -0.422                                          
## dstnc_ss_sc -0.065  0.178                                   
## daily_lngth -0.996  0.364  0.040                            
## foodlow     -0.426  0.076  0.032  0.409                     
## t_min_scald  0.366  0.039 -0.044 -0.389  0.144              
## precip_scld -0.016 -0.146  0.132  0.022  0.074 -0.153       
## arPlnttn:__  0.035  0.050 -0.649 -0.020 -0.059  0.049 -0.226
s2_dss_raw <- 
  ggplot(locs_day,
         aes(x = distance_ss,
             y = sum_S2,
             color = area)) +
  geom_point(
    aes(shape = area), 
    size = 3.2, 
    alpha = 0.3) +
  geom_smooth(
    aes(fill = area), method = "lm", 
              se = TRUE,
              linewidth = 2.8) +
  scale_color_manual(values = colores, 
                     name = "Study Site") +
  scale_fill_manual(values = colores2,
                    name = "Study Site") +
  scale_shape_manual(
    values = c("Plantations" = 16, "INP" = 17),
    name = "Study Site") +
  labs(
    x = "Distance Between Sleeping Sites (m)",
    y = "Daily Trajectories (m)") +
  guides(
    color = guide_legend(
      override.aes = list(
        shape = c(16, 17), 
        size = 5)),
    shape = "none") +
  theme_classic() +
  theme(
    text = element_text(family = 'Calibri', size = 34),
    legend.position = 'bottom',
    legend.box.margin = margin(t = 50, b = 10),
    legend.title = element_text(size = 34),
    legend.text = element_text(size = 30),
    legend.key.size = unit(1.5, "lines"),
    legend.spacing.x = unit(0.5, "cm"),
    legend.spacing.y = unit(0.3, "cm"),
    axis.title.x = element_text(size = 34, vjust = -3), 
    axis.title.y = element_text(size = 34, vjust = 3), 
    axis.text = element_text(size = 30),
    plot.margin = unit(c(1, 1, 1, 1), 'cm'))
s2_dss_raw

r_squared <- r.squaredGLMM(m21)

r_squared <- data.frame(
  Metric = c("Marginal R² (R²m)", "Conditional R² (R²c)"),
  Value = c(r_squared[1, 1], r_squared[1, 2]))

r_squared %>%
  flextable() %>%
  colformat_num(j = "Value", digits = 2) %>%
  set_caption("R² Values for best model - capuchins daily trajectories") %>%
  bold(part = "header") %>%
  autofit()
R² Values for best model - capuchins daily trajectories

Metric

Value

Marginal R² (R²m)

0.51

Conditional R² (R²c)

0.77

qqnorm(residuals(m21), cex = 1,
       cex.axis = 2.3,   
       cex.lab = 2.3,    
       cex.main = 2.3)

qqline(residuals(m21), col='darkred', cex=9)

plot(fitted(m21), residuals(m21),
     xlab = "Adjusted values", ylab = "Residuals",
     main = "Residuals vs. Adjusted values",
     cex.lab = 2.3,   
     cex.main = 2.3,   
     cex.axis = 2.3)

abline(h = 0, col = "red")

set.seed(123)


folds <- createFolds(locs_day$sum_S2, k = 10, list = TRUE)


errors <- numeric(length(folds))


for (i in seq_along(folds)) {
  train_data <- locs_day[-folds[[i]], ] 
  test_data <- locs_day[folds[[i]], ]  
  

  model <- lmer(sum_S2 ~ area * distance_ss_scaled + food + t_min_scaled + 
                  precip_scaled + (1 | grupo),
                data = train_data)  # Aquí usamos train_data
  

  predictions <- predict(model, newdata = test_data, allow.new.levels = TRUE)
  

  errors[i] <- mean((predictions - test_data$sum_S2)^2, na.rm = TRUE)
}

mean_cv_error <- mean(errors, na.rm = TRUE)
cat("RMSE cross-validation:", sqrt(mean_cv_error), "\n")
## RMSE cross-validation: 692.1672
model_no_random <- lm(sum_S2 ~ area * distance_ss_scaled + food + t_min_scaled + precip_scaled, 
                      data = locs_day)

errors_no_random <- numeric(length(folds))

for (i in seq_along(folds)) {
  train_data <- locs_day[-folds[[i]], ]
  test_data <- locs_day[folds[[i]], ]
  
  model <- lm(sum_S2 ~ area * distance_ss_scaled + food + t_min_scaled + precip_scaled, 
              data = train_data)
  
  predictions <- predict(model, newdata = test_data)
  
  errors_no_random[i] <- mean((predictions - test_data$sum_S2)^2, na.rm = TRUE)
}

mean_cv_error_no_random <- mean(errors_no_random, na.rm = TRUE)
cat("Mean RMSE without random effect:", sqrt(mean_cv_error_no_random), "\n")
## Mean RMSE without random effect: 692.2336

The linear mixed model indicates that capuchin group identity (random effect) accounts for some variation in daily travel distance, but the majority of the variance is explained by the fixed effects. The marginal R² (0.147) suggests that environmental and spatial predictors explain about 15% of the variation, while the conditional R² (0.149) shows that including the random effect adds little explanatory power. The RMSE values (692.17 vs. 692.23) are nearly identical, reinforcing that group-level variation does not substantially improve predictive accuracy. However, significant fixed effects—such as the interaction between study site and sleeping site distance, and greater travel distances in plantations—highlight the influence of habitat structure on capuchin movement patterns.

3.5 Plantaions groups - Bark stripping

locs_day_pine <- locs_day %>% filter(area == 'Plantations')

dano <- read_csv ('files/dano_nuevo_area.csv')

dano <- dano %>% 
  separate(mes, 
           into = c("mes", "day"), 
           sep = "-") 
dano <- dano %>% 
  rename(year= day)

dano <- dano  %>%
  mutate(year = paste0("20", year))

dano$grupo <- factor(dano$grupo, 
                     levels = c('Grupo 1', 'Grupo 2', 'Grupo 3', 'Grupo 4'), 
                     labels = c('Columba', 'Pandemia', 'Triangulo', 'Gualicho'))

monthly_dano <- dano %>%
  group_by(grupo, mes, year) %>%
  summarise(media_dano = mean(m2, na.rm = TRUE)) 

monthly_dano$year <- as.integer(monthly_dano$year)

monthly_dano <- monthly_dano %>% 
  rename(month = mes)

locs_day_pine$year <- as.integer(locs_day_pine$year)

locs_day_pine <- locs_day_pine %>%
  left_join(monthly_dano, by = c("grupo" = "grupo", "month" = "month", "year" = "year"))
ggplot(locs_day_pine, aes(x = media_dano, y = sum_S2)) +
  geom_jitter(
    size = 3, 
    alpha = 0.2, 
    color = 'black', 
    shape = 21, 
    fill = 'grey60',
    width = 0.08, 
    height = 0.1) +
  geom_smooth(
    method = "lm",
    linewidth = 1.8, 
    color = "midnightblue",
    fill = "lightblue",
    se = TRUE,
    alpha = 0.3) +
  labs(
    y = "Daily trajectories (m)",
    x = expression("Monthly bark-stripping per transect (m"^2*")")) +
  theme_bw() +
  guides(
    size = "none",
    alpha = "none",
    color = "none") +
  theme(
    text = element_text(family = "Calibri", size = 34),
    axis.text.x = element_text(color = "grey15", angle = 45, hjust = 0.5, vjust = 0.6),
    axis.text.y = element_text(color = "grey15"),
    axis.title.x = element_text(face = "bold", margin = margin(t = 20)),
    axis.title.y = element_text(face = "bold", vjust = 3),
    plot.margin = unit(c(1, 1, 1, 1), "cm"))

cor_test_result <- cor.test(locs_day_pine$media_dano, locs_day_pine$sum_S2, 
                            method = "pearson", use = "complete.obs")
print(cor_test_result)
## 
##  Pearson's product-moment correlation
## 
## data:  locs_day_pine$media_dano and locs_day_pine$sum_S2
## t = 1.8221, df = 329, p-value = 0.06934
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.007931638  0.205539895
## sample estimates:
##        cor 
## 0.09995422
# Extraer resultados relevantes
correlation_coefficient <- cor_test_result$estimate
p_value <- cor_test_result$p.value
n_obs <- cor_test_result$parameter + 2