To estimate capuchins monthly home ranges (95% isopleth) and core areas (50% isopleth) we fitted Ornstein–Uhlenbeck (OU) and Ornstein–Uhlenbeck Foraging (OUF) movement models. Datasets that better fit Independent Identically Distributed (IID), Brownian Motion (BM) were excluded because they did not accomplish the criteria described below. We selected the best-fitting model using the lowest Akaike Information Criterion (AIC; Akaike, 1974; Calabrese et al., 2016). After identifying the optimal movement model, we calculated Autocorrelated Kernel Density Estimates (AKDE), which account for positional and temporal autocorrelation. We specifically applied weighted AKDE to correct for uneven sampling, assigning higher weights to under-sampled intervals and lower weights to over-sampled intervals. This approach effectively handles datasets with missing or irregular data (Fleming et al., 2019; Noonan et al., 2019).

1 R Packages

library(sf)
library(ggplot2)
library(tidyverse)
library(gridExtra)
library(ctmm)
library(RColorBrewer)
library(patchwork)
library(lme4)
library(glmmTMB)
library(MuMIn)
library(broom.mixed)
library(flextable)
library(knitr)
library(lmerTest)
library(fitdistrplus)
library(cowplot)
library(lmtest)
library(caret)
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

Available on MoveBank (ID:)

hr <- read_rds('files/hr.RDS')

Summary

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

We have 21103 capuchins’ locations, 11348 for plantations landscape and 9755 for Iguazú National Park (INP). Plot for see data distribution

plantations_palette <- c("firebrick4","red","goldenrod2","darkorange","gold","sienna4")
inp_palette <- c("#023f09", "palegreen", "#2b9266", "green2", "olivedrab4")

hr <- hr %>%
  mutate(grupo = recode(grupo, "NI" = "Spot"))

hr <- hr %>%
  mutate(estacion2 = fct_recode(estacion2,
                             "Winter" = "invierno",
                             "Spring" = "primavera",
                             "Summer" = "verano",
                             "Autumn" = "otoño"))%>%
  rename(Group = grupo)

seasonal_counts <- hr %>%
  group_by(year, estacion2, Group, area) %>%
  summarise(count = n(), .groups = 'drop')

seasonal_counts$estacion2 <- factor(seasonal_counts$estacion2, 
                                    levels = c("Autumn", "Winter", "Spring", "Summer"))
y_limit <- 400

seasonal_counts <- seasonal_counts %>%
  mutate(high_value = ifelse(count > y_limit, "Above Limit", "Below Limit"))

# Plots
# INP

pni_sampling <- seasonal_counts %>%
  filter(area == 'PN Iguazu') %>%
  filter(!Group == 'Triangulo') %>% 
  ggplot(aes(x = estacion2, y = pmin(count, y_limit), fill = Group)) + 
  geom_col(position = position_dodge(width = 0.9), width = 1, color = "black", size = 1) +  
  geom_text(aes(label = ifelse(count > y_limit, paste0(" ", count), "")), 
            position = position_dodge(width = 0.8), 
            vjust = -0.5, size = 12, fontface = "bold") +  # Annotate bars > y_limit
  facet_wrap(~ year) +  
  labs(x = "Season", y = "Number of Locations") +
  scale_fill_manual(values = inp_palette) +  
  theme_minimal() +
  theme(
    text = element_text(family = "Calibri", size =60),
    axis.text.x = element_text(angle = 45, hjust = 1, face = "bold"),
    axis.text.y = element_text(),
    axis.title = element_text(, face = "bold"),
    strip.text = element_text(, face = "bold", color = "black", 
                              margin = margin(5, 5, 5, 5)),  
    strip.background = element_rect(fill = "grey74", color = "gray50"), 
    panel.grid.major = element_line(size = 0.3, color = "gray80"),
    panel.grid.minor = element_blank(),
    panel.spacing.y = unit(1, "lines"),
    legend.position = "right",
    legend.title = element_text(face = "bold"),
    legend.text = element_text(size = 45),
    legend.key.size = unit(0.9, "cm"),  
    legend.spacing.y = unit(1.2, "cm"), 
    plot.margin = margin(1, 1, 1, 1, "cm")) +
  ylim(0, y_limit + 50)

# Plantations

arauco_sampling <- seasonal_counts %>%
  filter(area == 'Arauco') %>% 
  ggplot(aes(x = estacion2, y = pmin(count, y_limit), fill = Group)) + 
  geom_col(position = position_dodge(width = 0.9), width = 1, color = "black", size = 1) +
  geom_bar(stat = "identity", position = "dodge") +  
  geom_text(aes(label = ifelse(count > y_limit, paste0(" ", count), "")), 
            position = position_dodge(width = 1), size = 12, angle =45, vjust = -0.05, hjust = -0.05) +  
  facet_wrap(~ year) +  
  labs(x = "Season", y = "Number of Locations") +
  scale_fill_manual(values = plantations_palette) +  
  theme_minimal() +
  theme(
    text = element_text(family = "Calibri", size = 60),
    axis.text.x = element_text(angle = 45, hjust = 1,  face = "bold"),
    axis.text.y = element_text(),
    axis.title.x = element_text(size = 0),
    axis.title = element_text(face = "bold"),
    strip.text = element_text( face = "bold", color = "black", 
                              margin = margin(5, 5, 5, 5)),  
    strip.background = element_rect(fill = "grey74", color = "gray50"), 
    panel.grid.major = element_line(size = 0.3, color = "gray80"),
    panel.grid.minor = element_blank(),
    panel.spacing.y = unit(1, "lines"),
    legend.position = "right",
    legend.title = element_text(face = "bold"),
    legend.text = element_text(size = 45),
    legend.key.size = unit(0.9, "cm"),  
    legend.spacing.y = unit(1.2, "cm"), 
    plot.margin = margin(1, 1, 1, 1, "cm")) +
  ylim(0, y_limit + 50)

arauco_sampling <- arauco_sampling + 
  labs(title = "Plantations") +
  theme(plot.title = element_text(hjust = 0, face = "bold"))

pni_sampling <- pni_sampling + 
  labs(title = "INP") +
  theme(plot.title = element_text( hjust = 0, face = "bold"))

# Panel plot
fs1_plot <- (arauco_sampling / pni_sampling) + 
  plot_layout(heights = c(0.7, 1.6)) + 
  plot_annotation(
    theme = theme(
      plot.title = element_text( face = "bold", hjust = 0.5)))
fs1_plot
*Figure S1. Locations of Black capuchins (_Sapajus nigritus_) groups in pine plantations and Iguazú National Park (INP) in northeastern Argentina.*

Figure S1. Locations of Black capuchins (Sapajus nigritus) groups in pine plantations and Iguazú National Park (INP) in northeastern Argentina. 

ggsave("fs1_plot_highres.png",
       plot = fs1_plot,
       width = 4500,
       height = 3900,
       units = "px",
       dpi = 300)

3 Calculating AKDEs

The data in the following table was obtained after running a function I created: ‘sim_akde.R’.
For more details on the function, check the link.

This function generates a summary table of the 95% and 50% AKDEs for each independent dataset; in this case, the monthly locations for each group (AKDEs from different months).

Datasets with small effective size (<30) were discarded. We discarted AKDEs with no residence behavior based on variograms plotted using ‘ctmm’ package.
I do not think it is appropriate to compare areas estimated under an autocorrelated movement model (which adjusts the bandwidth based on point correlation; OU, OUF) and non-autocorrelated models.

Additionally, sim_akde.R generates a shapefile for each AKDE, which we will use later to estimate proportions of each habitat type within plantations landscapes (pine stands and native forest fragments).

3.0.1 Rarerifying data

rarefy_gps <- function(data) {
  data <- data %>%
    mutate(
      date = as.Date(timestamp),
      rounded_time = timestamp - (as.numeric(timestamp) %% (90 * 60))
    ) %>%
    group_by(individual.local.identifier, date) %>%
    filter(abs(difftime(timestamp, rounded_time, units = "mins")) <= 15) %>%
    slice_min(order_by = abs(difftime(timestamp, rounded_time, units = "mins")), n = 8, with_ties = FALSE) %>%
    ungroup() %>%
    arrange(individual.local.identifier, timestamp)
  
  return(data)
}

hr_stand <- rarefy_gps(hr)
sim_akde <- function(data, output_dir = "ctmm_fits/monthly/") {
  
  shapes_dir <- file.path(output_dir, "shapes")
  
  if(!dir.exists(output_dir)) dir.create(output_dir, recursive = TRUE)
  if(!dir.exists(shapes_dir)) dir.create(shapes_dir, recursive = TRUE)
  
  results <- data.frame(
    individual.local.identifier = character(),
    best_fit = character(),
    akde_95_area = numeric(),
    akde_50_area = numeric(),
    n_locs = integer(),
    stringsAsFactors = FALSE
  )
  
  unique_ids <- unique(data$individual.local.identifier)
  
  for (id in unique_ids) {
    
    telemetry_data <- data %>% filter(individual.local.identifier == id)
    
    telemetry_object <- as.telemetry(telemetry_data, timezone = "America/Argentina/Buenos_Aires")
    
    guess <- ctmm.guess(telemetry_object, CTMM = ctmm(error = 20), interactive = FALSE)
    
    fits <- ctmm.select(telemetry_object, CTMM = guess, verbose = TRUE)
    
    saveRDS(fits, file = file.path(output_dir, paste0(id, "_fits.RDS")))
    
    akde_result <- akde(data = telemetry_object, CTMM = fits[[1]], weights = TRUE)
    
    n_eff <- summary(akde_result)$DOF["n_eff"]
    if (n_eff > 30) {
      akde_sf_95 <- as.sf(akde_result, level.UD = 0.95, level = 0.95) %>% st_transform(32721)
      akde_sf_50 <- as.sf(akde_result, level.UD = 0.5, level = 0.5) %>% st_transform(32721)
      
      akde_95 <- akde_sf_95 %>% filter(name == paste0(id, " 95% est"))
      akde_50 <- akde_sf_50 %>% filter(name == paste0(id, " 50% est"))
      
      st_write(akde_95, file.path(shapes_dir, paste0(id, "_akde_95.shp")), delete_layer = TRUE)
      st_write(akde_50, file.path(shapes_dir, paste0(id, "_akde_50.shp")), delete_layer = TRUE)
      
      best_fit_name <- if (is.list(fits) && length(fits) > 0) {
        summary(fits[[1]])$name
      } else {
        NA
      }
      
      results <- rbind(results, data.frame(
        individual.local.identifier = id,
        best_fit = best_fit_name,
        akde_95_area = as.numeric(st_area(akde_95)),
        akde_50_area = as.numeric(st_area(akde_50)),
        n_locs = nrow(telemetry_data)
      ))
    }
  }
  
  return(results)
}

Apply the function

hr_month <- sim_akde(hr_stand)
hr_month <- hr_month %>% 
  mutate(across(c(area,
                  best_fit),
                as.factor))

hr_month <- hr_month %>%
  separate(individual.local.identifier, 
           into = c("grupo", "year", "mes"), 
           sep = "_") %>% 
  mutate(ratio = akde_95_area/akde_50_area)
summary_table <- hr_month %>%
  group_by(area) %>%
  summarise(
    mean_akde95 = mean(akde_95_area, na.rm = TRUE),
    sd_akde95 = sd(akde_95_area, na.rm = TRUE),
    mean_akde50 = mean(akde_50_area, na.rm = TRUE),
    sd_akde50 = sd(akde_50_area, na.rm = TRUE),
    mean_ratio = mean(akde_95_area / akde_50_area, na.rm = TRUE),
    sd_ratio = sd(akde_95_area / akde_50_area, na.rm = TRUE)) %>%
  mutate(area = recode(area, "arauco" = "Plantations", "pni" = "INP"))

ft <- flextable(summary_table) %>%
  set_header_labels(
    mean_akde95 = "Mean AKDE 95% (ha)",
    sd_akde95 = "SD AKDE 95%",
    mean_akde50 = "Mean AKDE 50% (ha)",
    sd_akde50 = "SD AKDE 50%",
    mean_ratio = "Mean Ratio (95% / 50%)",
    sd_ratio = "SD Ratio"
  ) %>%
  theme_tron_legacy() %>%
  colformat_num(j = c("mean_akde95", "sd_akde95", "mean_akde50", "sd_akde50", "mean_ratio", "sd_ratio"), digits = 2, x = .) %>%
  autofit() %>%
  set_caption("Summary of Home Range and Core Areas")

ft
Summary of Home Range and Core Areas

area

Mean AKDE 95% (ha)

SD AKDE 95%

Mean AKDE 50% (ha)

SD AKDE 50%

Mean Ratio (95% / 50%)

SD Ratio

Plantations

392.0325

71.94982

72.2628

27.10073

6.074641

2.1424860

INP

134.1412

45.18878

41.7681

18.17756

3.408943

0.8244173

4 Covariates

4.1 Landscape ShapeFiles

The sim_akde function saved multiple Shapefiles corresponding to each home range (AKDE95) and core area (AKDE50). These shapefiles are not only useful for creating maps but also for calculating important metrics, such as the proportion of Native Forest without management (BNSM; abreviation in spanish) or pine plantations.

There are many shapefiles (2 per home range, totaling 190).

I read all the Shapefiles from the ‘monthly’ folder (where sim_akde saved them) based on their file extension. I kept the full.names option because sim_akde saved the files preserving the group_year_month format, which is useful to later split them into a dataframe and match them with the summary table of the AKDEs.

Applying st_read to all elements in the list was problematic, since it generated too many automatic messages in the HTML output, including CRS information for each of the 190 shapefiles, bounding boxes, etc. To avoid this, I did the following:

files_shp <-
  list.files(path = "/shapes",
             pattern = "\\.shp$",
             full.names = TRUE)

lista_shapefiles <- lapply(files_shp, st_read)

Need to add an aditional shape of other timber company

BNSM_shapefile <- st_read("files/BNSM.shp")
## Reading layer `bnsm' from data source 
##   `C:\Users\Usuario\Desktop\Doctorado\Cap2\paper\codes\paper_codes\files\bnsm.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 217 features and 21 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -54.53957 ymin: -25.93148 xmax: -54.26311 ymax: -25.74938
## Geodetic CRS:  WGS 84

To store all the results generated during the process, I first create an empty list.

The first line extracts the CRS used by Arauco for their shapefile (in case they used an uncommon one).

lapply iterates over all the shapefiles (190 in total) and reprojects the coordinates according to the CRS of the BNSM layer using st_transform(). This ensures that all shapefiles in lista_shapefiles are in the same CRS.

  • The for loop iterates through each shapefile in lista_shapefiles_transformada.
  • In each iteration, it extracts the filename using basename(), which returns only the file name without the full path (useful for identifying which file is being processed).
  • Then, it calculates the geometric intersection between the transformed shapefile and BNSM_shapefile using st_intersection(). This returns the geometry corresponding to the overlapping area between the two shapefiles.
  • Next, it calculates the total area of the transformed shapefile using st_area(), which returns the areas of all geometries in the shapefile, and then sums these areas using sum().
  • It also calculates the area of the intersection (i.e., the overlapping region between the current shapefile and BNSM_shapefile).
  • Finally, it computes the proportion by dividing the intersection area by the total area and stores all the results in a data frame.
resultados <- list()
crs_bsnm <- st_crs(BNSM_shapefile)

lista_shapefiles_transformada <- lapply(lista_shapefiles, function(shapefile) {
  st_transform(shapefile, crs_bsnm)
})


for (i in seq_along(lista_shapefiles_transformada)) {

  nombre_shape <- basename(archivos_shp[i])
  

  interseccion <- st_intersection(lista_shapefiles_transformada[[i]], BNSM_shapefile)
  

  area_shape <- sum(st_area(lista_shapefiles_transformada[[i]]))
  area_interseccion <- sum(st_area(interseccion))
  
 
  proporcion <- as.numeric(area_interseccion / area_shape)
  
  
  resultados[[i]] <- data.frame(
    Shape = nombre_shape,
    Proporcion_BNSM = proporcion)
}

df_resultados <- do.call(rbind, resultados)


df_resultados <- df_resultados %>%
  separate(Shape, 
           into = c('Group', 'year', 'mes', 'akde', 'isopleth'), 
           sep = '_')

df_resultados_wide <- df_resultados %>%
  mutate(isopleth = gsub("\\.shp", "", isopleth)) %>%  
  pivot_wider(names_from = isopleth, 
              values_from = Proporcion_BNSM, 
              names_prefix = "prop_bnsm_")

hr_month <- hr_month %>%
  left_join(df_resultados_wide, by = c("grupo", "year", "mes"))

4.1.1 Checking spatial information

BNSM_shapefile <- st_transform(BNSM_shapefile, crs = 4326)  
lista_shapefiles[[3]] <- st_transform(lista_shapefiles[[3]], crs = 4326)

# CRS 

lista_shapefiles[[3]] <- st_transform(lista_shapefiles[[3]], st_crs(BNSM_shapefile))

# bbox for plot limits, we are using shape '3' randomnly as a test

bbox_shape <- st_bbox(lista_shapefiles[[3]])

# adjust a little bit 
x_margin <- (bbox_shape["xmax"] - bbox_shape["xmin"]) * 0.4
y_margin <- (bbox_shape["ymax"] - bbox_shape["ymin"]) * 0.4
ggplot() +
  geom_sf(
    data = BNSM_shapefile,
    fill = "olivedrab3",
    color = "forestgreen") +
  geom_sf(
    data = lista_shapefiles[[3]],
    fill = adjustcolor("goldenrod1", alpha.f = 0.15),
    color = "gold",
    lwd = 1.5) +
  coord_sf(
    xlim = c(bbox_shape["xmin"] - x_margin, bbox_shape["xmax"] + x_margin),
    ylim = c(bbox_shape["ymin"] - y_margin, bbox_shape["ymax"] + y_margin),
    expand = FALSE) +
  theme_minimal() +
  labs() +
  theme(
    text = element_text(
      family = "Calibri",
      face = "bold",
      size = 34),
    axis.text.x = element_text(
      color = "black", 
      angle =45,
      vjust =0.5,
      face = "bold"),
    axis.text.y = element_text(
      color = "black",
      face = "bold"),
    panel.background = element_rect(fill = "saddlebrown"),
    plot.background = element_rect(fill = "white", color = NA),
    panel.grid.major = element_line(color = "saddlebrown"),
    panel.grid.minor = element_line(color = "saddlebrown"))
Capuchin home ranges in pine plantations. Plantations stands are in brown and native forest fragments are in green.

Capuchin home ranges in pine plantations. Plantations stands are in brown and native forest fragments are in green.

Lets see from which group was that AKDE, and see if it make sense

lista_shapefiles[[3]]$name
## [1] "Columba_2022_abr 95% est"
ej <- hr_month %>% 
  filter(grupo == "Columba", year == "2022", mes == "abr") %>% 
  pull(prop_bnsm_95)

paste("prop_bnsm_95:", ej)
## [1] "prop_bnsm_95: 0.189329491"

4.2 Bark-stripping

Over a 2-year sampling period (2021–2023), we monthly measured intensity of bark-stripping on 2000 pine trees distributed across 40 transects, randomly placed within the home ranges of four black capuchin groups. As a proxy of pine phloem consumption, we measured the total bark-stripped area per pine, calculated as the lateral surface area of a truncated cylinder (Kershaw et al. 2016), using the recorded start and end heights of the debarked area and a taper parameter (k), derived from preliminary measurements on 24 pines to adjust for tree growth patterns specific to the study area. \[ A = \pi \times \frac{D_i}{2} \left( \frac{1}{1+ e^{k (h_{\text{lower}}- h_i ) }} + \frac{1}{1+ e^{k (h_{\text{upper}}- h_i ) }} \right) \times \sqrt{ h^2 + \frac{D_i}{2} \left( \frac{1}{1+ e^{k (h_{\text{lower}}- h_i ) }} + \frac{1}{1+ e^{k (h_{\text{upper}}- h_i ) }} \right) } \] Where \(A\) is the total lateral surface area debarked,
\(D_i\) is the diameter at breast height (DBH; 1.3 m above the ground),
\(k\) is the taper parameter that adjusts the rate of diameter reduction with height,
\(h_{\text{lower}}\) and \(h_{\text{upper}}\) are the heights where the bark-stripping starts and ends, respectively,
and \(h_i\) is the fixed reference height (1.3 m).

We will add the bark-stripping monitoring to the AKDE dataset

dano <- readRDS('files/dano_paper_final.rds')

Plotting the monthly bark-stripping pattern through the two-year cycle

ggplot(dano, aes(x = month, y = m2)) + 
 
  geom_jitter(width = 0.1, alpha = 0.04, shape = 16, size = 5) +
  
  geom_smooth(
    data = dano,
    aes(x = month, y = m2, group = edad, color = edad, linetype = edad), 
    method = "loess",
    se = FALSE,
    size = 2.2
  ) +
  
  geom_smooth(
    data = dano, 
    aes(group = 1, color = "Average", linetype = "Average"), 
    method = "loess", 
    se = TRUE, 
    size = 2.6
  ) +
  
  labs(
    x = "Month", 
    y = bquote("Monthly bark-stripping" ~ (m^2)) ) + 
  
  scale_x_discrete(labels = c("Jan", "Feb", "Mar", "Apr", "May", "Jun", 
                              "Jul", "Aug", "Sep", "Oct", "Nov", "Dec")) + 
  
  scale_y_continuous(limits = c(0, 3)) +

  scale_color_manual(values = c(
    "Average" = "#401919",  
    "mature pine stands" = "#d07463",
    "medium-aged pine stands" = "#cd9253",
    "young pine stands" = "#ced063")) +
  scale_linetype_manual(values = c(
    "Average" = "solid",
    "mature pine stands" = "dashed",
    "medium-aged pine stands" = "dashed",
    "young pine stands" = "dashed")) +
  theme_classic() + 
  theme(
    text = element_text(
      family = "Calibri",
      size = 34),
    axis.title.x = element_text(
      margin = margin(t = 17),
      color = "grey15"),
    axis.text.x = element_text(angle = 45, 
                               color = "grey15", 
                               vjust = 0.6,
                               hjust = 0.5),
    axis.title.y = element_text(margin = margin(r = 5),
                                color = "grey15"),
    axis.text.y = element_text(color = 'grey15'),
    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"),
    legend.position = "bottom", 
    legend.title = element_text(size = 0),
    legend.text = element_text(size = 28)) +
  guides(
    color = guide_legend(
      nrow = 2,               
      byrow = TRUE,          
      override.aes = list(    
        linetype = c( "solid", "dashed", "dashed", "dashed"),  
        size = c(1.8, 1.8, 1.8, 2.6))),
    linetype = guide_legend(
      nrow = 2,
      byrow = TRUE))
*Figure B.2.: Monthly bark-stripped area recorded during a two-year survey (2021–2023) of ~2,000 pine trees categorized by age classes.*

Figure B.2.: Monthly bark-stripped area recorded during a two-year survey (2021–2023) of ~2,000 pine trees categorized by age classes.

4.3 Fruit Ripening

For capuchins living in plantations and Iguazú National Park (INP), we included fruit ripening season as a covariate, as it has been reported to influence space use patterns (Di Bitetti 2001a) . Since our data spans over a decade, we established a broad categorization based on two independent data sources. To define this classification, we used phenology data from Di Bitetti (2001a), collected in INP, and 2021–2023 phenological surveys conducted in San Jorge Reserve, a private protected area located between the plantation study site and INP (-25.766703, -54.309295). Phenological surveys conducted between June 2021 and October 2023 consisted of transect-based monitoring, in which we bimonthly recorded the phenological state of 416 individuals from 52 tree and liana species known to be consumed by capuchins (Di Bitetti 2001b). During each survey, we classified fruit ripening into four levels: 0 (no ripe fruits), 1 (0–10% ripening relative to its maximum), 2 (10–30%), and 3 (>30%), following Agostini et al. (2010). Each species was further categorized according to its importance in the capuchin diet (Di Bitetti 2001b) as Low, Intermediate, or High. Based on this information (Figure B.1), we classified May–August as low fruit ripening months and September–April as high fruit ripening months.

smoothed_feno <- readRDS('files/smoothed_feno.rds')

color_palette <- c("High" = "#481567",  
                   "Intermediate" = "#21908C",  
                   "Low" = "#35B779")  


smoothed_feno <- smoothed_feno %>%
  mutate(imp_label = factor(imp, levels = c(1, 2, 3), labels = c("Low", "Intermediate", "High")))

phenology_paper <-
  ggplot(smoothed_feno, aes(x = month, y = smoothed_mean, color = imp_label, group = imp_label)) +
  geom_ribbon(aes(ymin = smoothed_mean - smoothed_sd, 
                  ymax = smoothed_mean + smoothed_sd, 
                  fill = imp_label), 
              alpha = 0.15, color = NA) + 
  geom_smooth(method = "loess", se = FALSE, size = 1.9, span = 0.6) +  
  theme_classic() +
  scale_x_discrete(limits = month.name) +  
  scale_color_manual(values = color_palette, name = "Importance in Diet") +  
  scale_fill_manual(values = color_palette, name = "Importance in Diet") +   
  labs(y = "Mean Fruit Rippening", x = "Month") +  
  theme(
    text = element_text(family = "Calibri", size = 34),
    axis.text.x = element_text(angle = 45, hjust = 1),
    axis.title.x = element_text(vjust = 0.5),
    axis.title.y = element_text(vjust = 2.5),
    legend.position = "top", 
    legend.text = element_text(size = 30),  
    legend.title = element_text(face = "bold", size = 32), 
    plot.margin = margin(1, 1, 1, 1, unit = "cm")  
  ) +
  guides(color = guide_legend(title = "Importance in Diet"),
         fill = guide_legend(title = "Importance in Diet"))
  
phenology_paper
*Figure B.1.: Seasonal trends in fruit ripening levels based on phenological surveys conducted between June 2021 and October 2023 in northeastern Argentina. The solid lines represent the mean fruit ripening index for tree and liana species consumed by Black Capuchins Monkeys (Sapajus nigritus), categorized by their dietary importance: High (purple), Intermediate (teal), and Low (green). The shaded areas indicate the standard deviation for each category. Fruit ripening was classified into four levels: 0 (no ripe fruits), 1 (0–10% ripening), 2 (10–30%), and 3 (>30%). Based on these patterns, we classified May–August as low fruit ripening months and September–April as high fruit ripening months.*

Figure B.1.: Seasonal trends in fruit ripening levels based on phenological surveys conducted between June 2021 and October 2023 in northeastern Argentina. The solid lines represent the mean fruit ripening index for tree and liana species consumed by Black Capuchins Monkeys (Sapajus nigritus), categorized by their dietary importance: High (purple), Intermediate (teal), and Low (green). The shaded areas indicate the standard deviation for each category. Fruit ripening was classified into four levels: 0 (no ripe fruits), 1 (0–10% ripening), 2 (10–30%), and 3 (>30%). Based on these patterns, we classified May–August as low fruit ripening months and September–April as high fruit ripening months.

5 AKDEs GLMM

5.1 AKDE 95% (Home ranges)

We preserved AKDe that overcome the criteria of having effective size over 30 (indepedent locations) and ranging behavior (based on variograms). Summarizing data, AKDE were

A summary of AKDEs before the statistical analysis

best_fit_summary <- hr_month %>%
  count(best_fit) %>%
  mutate(percentage = (n / sum(n)) * 100)

hr_count_per_month <- hr_month %>%
  count(month, area)
colores <- c('orange3', 'palegreen4')

ggplot(hr_count_per_month, aes(x = month, y = n, fill = area)) +
  geom_bar(stat = "identity", position = "dodge", color = "black", alpha = 0.9) +
  scale_fill_manual(values = colores) +
  labs(title = "AKDEs per Month by Area",
       x = "Month",
       y = "Number of Home Ranges",
       fill = "Study Area") +
  theme_minimal() +
  theme(
    text = element_text(family = "Calibri", size =34) ,
    axis.text.x = element_text(angle = 45, hjust = 1))

hr_month <- hr_month %>%
  mutate(month = recode(mes,
                        "ene" = "January", "feb" = "February", "mar" = "March",
                        "abr" = "April", "may" = "May", "jun" = "June",
                        "jul" = "July", "ago" = "August", "sep" = "September",
                        "oct" = "October", "nov" = "November", "dic" = "December"))
normal=fitdist(hr_month$akde_95_area,"norm")

lognormal=fitdist(hr_month$akde_95_area,"lnorm")

gamma=fitdist(hr_month$akde_95_area,"gamma")


cdf.akde = cdfcomp(
  list(lognormal, normal, gamma),
  main = "",
  legendtext = c("lognormal", "normal", "gama"),
  fitcol = c("darkred", "orange3", "goldenrod1"),
  plotstyle = "ggplot") +
  geom_line(size = 1.6) +
  theme(
    text=element_text(size = 34),
    axis.title = element_text(),
    axis.text = element_text(),
    axis.text.x = element_text(angle = 45,
                               vjust = 0.8),
    legend.position = c(1.6, 0.9),
    legend.text = element_text())

qq.akde = qqcomp(
  list(lognormal, normal, gamma),
  main = "",
  fitcol = c("darkred", "orange3", "goldenrod1"),
  plotstyle     = "ggplot") +
  geom_line(size = 1.6) +
  theme_bw() +
  theme(
    text=element_text(size = 34),
    axis.text.x = element_text(angle = 45,
                               vjust = 0.8),
    legend.position = "none")

plot_grid(qq.akde, cdf.akde, NULL, ncol = 3, rel_widths = c(1.5, 1.5, 0.75))

gofstat(list(normal, lognormal, gamma),
        fitnames = c("normal","Lognormal","gamma"))[12]
## $aic
##    normal Lognormal     gamma 
##  1213.263  1193.296  1194.366

We will go with lognormal distribution on respose variable (akde95)

hr_month$log_akde_95_area <- log(hr_month$akde_95_area)

m0 <- lmer(log_akde_95_area ~ 1 + (1 | grupo), data = hr_month)
m1 <- lmer(log_akde_95_area ~ area + (1 | grupo), data = hr_month)
m2 <- lmer(log_akde_95_area ~ area + productividad + (1 | grupo), data = hr_month)
m3 <- lmer(log_akde_95_area ~ area + productividad + n_ind + (1 | grupo), data = hr_month)
m4 <- lmer(log_akde_95_area ~ area + n_ind + (1 | grupo), data = hr_month)
m5 <- lmer(log_akde_95_area ~ area * productividad + n_ind + (1 | grupo), data = hr_month)
m6 <- lmer(log_akde_95_area ~ area * n_ind + productividad + (1 | grupo), data = hr_month)
m7 <- lmer(log_akde_95_area ~ area * n_ind * productividad + (1 | grupo), data = hr_month)
AIC(m0,m1,m2,m3, m4, m5, m6, m7)
glmm_summary <- tidy(m3, effects = "fixed") %>%
  mutate(p.value = format.pval(p.value, digits = 3, eps = 0.001)) %>%
  rename(
    Variable = term,              
    Estimate = estimate,          
    `Standard_Error` = std.error,  
    `P_Value` = p.value) %>%

  mutate(Variable = recode(Variable,
    "area" = "Study site (Plantations vs. INP)", 
    "productividad" = "Fruit Ripening",
    "n_ind" = "Group size"))


glmm_results %>%
  flextable() %>%
  theme_tron_legacy() %>%
  colformat_num(j = c("Estimate", "Standard_Error", "P_Value"), digits = 3) %>%  
  bold(part = "header") %>%
  autofit() %>%
  set_caption("Best fitted GLMM, Home ranges (AKDE 95%)")
qqnorm(residuals(m3), cex = 1,
       cex.axis = 2.3,   
       cex.lab = 2.3,    
       cex.main = 2.3)

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

shapiro.test(residuals(m3))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(m3)
## W = 0.98962, p-value = 0.6684
plot(fitted(m3), residuals(m3),
   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")

residuos <- residuals(m3)

valores_ajustados <- fitted(m3)

modelo_aux <- lm(residuos^2 ~ valores_ajustados)

bptest(modelo_aux)
## 
##  studentized Breusch-Pagan test
## 
## data:  modelo_aux
## BP = 1.2554, df = 1, p-value = 0.2625
r_squared <- r.squaredGLMM(m3)
r2m <- r_squared[1, 1] 
r2c <- r_squared[1, 2] 


set.seed(123)
folds <- createFolds(hr_month2$log_akde_95_area, k = 10, list = TRUE)
errors <- numeric(length(folds))

for (i in seq_along(folds)) {
  train_data <- hr_month2[-folds[[i]], ]  
  test_data <- hr_month2[folds[[i]], ]   
  
  model <- lmer(log_akde_95_area ~ area + productividad + n_ind + (1 | grupo), 
                data = train_data)
  
  predictions <- predict(model, newdata = test_data, allow.new.levels = TRUE)
  
  errors[i] <- mean((predictions - test_data$log_akde_95_area)^2)
}

mean_cv_error <- mean(errors)
rmse_log <- sqrt(mean_cv_error)  
typical_error_percentage <- (exp(rmse_log) - 1) * 100  


model_performance <- data.frame(
  Metric = c("Marginal R² (R²m)", "Conditional R² (R²c)", "RMSE (log scale)", "Typical Error (%)"),
  Value = c(r2m, r2c, rmse_log, typical_error_percentage)
)


model_performance %>%
  flextable() %>%
  theme_tron_legacy() %>%  
  colformat_num(j = "Value", digits = 3) %>%
  bold(part = "header") %>%
  autofit() %>%
  set_caption("Model Performance Summary: R² and RMSE")
Model Performance Summary: R² and RMSE

Metric

Value

Marginal R² (R²m)

0.7268485

Conditional R² (R²c)

0.9341781

RMSE (log scale)

0.1966613

Typical Error (%)

21.7331704

The selected model explained a substantial proportion of the variance in home range size, with fixed effects accounting for 78.85% (R²m = 0.789) of the variance and the total model explaining 90.51% (R²c = 0.905).

Cross-validation supported the robustness of the model, yielding an average root mean square error (RMSE) of 0.1966 in the log-transformed scale. When back-transformed to the original scale, this corresponds to an approximate 2.17 × 10¹% typical error.

colores <- c('orange3', 'palegreen4')


ggplot(hr_month, aes(x = area, y = akde_95_area, fill = area)) +
    facet_wrap(~ productividad) +
    geom_boxplot(outlier.shape = NA, lwd = 1.1, color = 'black', alpha = 0.8) +
    scale_fill_manual(values = colores) +
    geom_jitter(
      aes(size = n_ind),
      alpha = 0.4,
      shape = 21,
      color = 'black',
      width = 0.4,
      height = 5) + 
    scale_y_continuous(breaks = seq(60, 520, by = 40),
                       limits = c(60, 520)) +
    theme_bw() +
    guides(fill = "none",
           color = "none",
           alpha = 'none') +
    labs(
      y = "AKDE 95% (ha)",
      title = "Fruit ripening",
      fill = "Study site",
      size = "Group Size") +
    theme(
      text = element_text(
        family = "Calibri",
        size = 34),
      axis.text.x = element_text(
        size = 0,
        angle = 55,
        color = "grey18",
        vjust = 0.6,
        hjust = 0.5),
      axis.title.x = element_text(size = 0),
      axis.text.y = element_text(color = 'grey15'),
      axis.title.y = element_text(vjust = 3),
      legend.box = element_blank(),
      plot.margin = unit(c(1, 1, 1, 1), "cm"))
*Figure 2: Predicted monthly Auto-correlated Kernel Density Estimator (AKDE) 95% (Home Ranges) of black capuchins (Sapajus nigritus) in pine plantations (orange) and Iguazú National Park (INP, green) across low and high ripe fruit availability seasons. Home ranges were similar in both study sites across fruit ripening seasons*

Figure 2: Predicted monthly Auto-correlated Kernel Density Estimator (AKDE) 95% (Home Ranges) of black capuchins (Sapajus nigritus) in pine plantations (orange) and Iguazú National Park (INP, green) across low and high ripe fruit availability seasons. Home ranges were similar in both study sites across fruit ripening seasons

5.2 AKDE 50% (Core areas)

normal=fitdist(hr_month$akde_50_area,"norm")

lognormal=fitdist(hr_month$akde_50_area,"lnorm")

gamma=fitdist(hr_month$akde_50_area,"gamma")


cdf.akde = cdfcomp(
  list(lognormal, normal, gamma),
  main = "",
  legendtext = c("lognormal", "normal", "gama"),
  fitcol = c("darkred", "orange3", "goldenrod1"),
  plotstyle = "ggplot") +
  geom_line(size = 1.6) +
  theme(
    text =element_text(size =34),
    axis.text.x = element_text(angle = 45,
                               vjust = 0.8),
    legend.position = c(1.6, 0.9),
    legend.text = element_text(size = 20))

qq.akde = qqcomp(
  list(lognormal, normal, gamma),
  main = "",
  fitcol = c("darkred", "orange3", "goldenrod1"),
  plotstyle     = "ggplot") +
  geom_line(size = 1.6) +
  theme_bw() +
  theme(
    text =element_text(size =34),
    axis.text.x = element_text(angle = 45,
                               vjust = 0.8),
    legend.position = "none")

plot_grid(qq.akde, cdf.akde, NULL, ncol = 3, rel_widths = c(1.5, 1.5, 0.75))

gofstat(list(normal, lognormal, gamma),
        fitnames = c("normal","Lognormal","gamma"))[12]
## $aic
##    normal Lognormal     gamma 
##  899.5726  877.3327  879.4907
hr_month$log_akde_50_area <- log(hr_month$akde_50_area)

m0_50 <- lmer(log_akde_50_area ~ 1 + (1 | grupo), data = hr_month)
m1_50 <- lmer(log_akde_50_area ~ area + (1 | grupo), data = hr_month)
m2_50 <- lmer(log_akde_50_area ~ area + productividad + (1 | grupo), data = hr_month)
m3_50 <- lmer(log_akde_50_area ~ area + productividad + n_ind + (1 | grupo), data = hr_month)
m4_50 <- lmer(log_akde_50_area ~ area + n_ind + (1 | grupo), data = hr_month)
m5_50 <- lmer(log_akde_50_area ~ area * productividad + n_ind + (1 | grupo), data = hr_month)
m6_50 <- lmer(log_akde_50_area ~ area * n_ind + productividad + (1 | grupo), data = hr_month)
m7_50 <- lmer(log_akde_50_area ~ area * n_ind * productividad + (1 | grupo), data = hr_month)
AIC(m0_50,m1_50,m2_50,m3_50,m4_50, m5_50, m6_50, m7_50)
glmm_summary <- tidy(m3_50, effects = "fixed") %>%
  mutate(p.value = format.pval(p.value, digits = 3, eps = 0.001)) %>%
  rename(
    Variable = term,              
    Estimate = estimate,          
    `Standard_Error` = std.error,  
    `P_Value` = p.value) %>%

  mutate(Variable = recode(Variable,
    "area" = "Study site (Plantations vs. INP)", 
    "productividad" = "Fruit Ripening",
    "n_ind" = "Group size"))


glmm_results %>%
  flextable() %>%
  theme_tron_legacy() %>%
  colformat_num(j = c("Estimate", "Standard_Error", "P_Value"), digits = 3) %>%  
  bold(part = "header") %>%
  autofit() %>%
  set_caption("Best fitted GLMM, Core areas (AKDE 50%)")
glmm_results <- data.frame(
  Variable = c("Study site (Plantations vs. INP)", "Group size", "Fruit Ripening"),
  Estimate = c(-0.48, 0.067, 0.193),
  `Standard_Error` = c(0.34, 0.013, -0.193),  # Rename to avoid spaces
  `P_Value` = c(0.103, 0.0604, 0.0345)
)

glmm_results %>%
  flextable() %>%
  theme_tron_legacy() %>%
  colformat_num(j = c("Estimate", "Standard_Error", "P_Value"), digits = 3) %>%  # Match exact column names
  bold(part = "header") %>%
  autofit() %>%
  set_caption("Best fitted GLMM, Core areas (AKDE 50%)")
Best fitted GLMM, Core areas (AKDE 50%)

Variable

Estimate

Standard_Error

P_Value

Study site (Plantations vs. INP)

-0.480

0.340

0.1030

Group size

0.067

0.013

0.0604

Fruit Ripening

0.193

-0.193

0.0345

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

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

plot(fitted(m3_50), residuals(m3_50),
     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")

residuos <- residuals(m3_50)

valores_ajustados <- fitted(m3_50)

modelo_aux <- lm(residuos^2 ~ valores_ajustados)

bptest(modelo_aux)
## 
##  studentized Breusch-Pagan test
## 
## data:  modelo_aux
## BP = 2.1801, df = 1, p-value = 0.1398
colores <- c('orange3', 'palegreen4')

hr_month <- hr_month %>%
  mutate(productividad = factor(productividad, levels = c("Baja", "Alta"), labels = c("Low", "High")))

ggplot(hr_month, aes(x = area, y = akde_50_area, fill = area)) +
    facet_wrap(~ productividad) +
    geom_boxplot(outlier.shape = NA, lwd = 1.1, color = 'black', alpha = 0.8) +
    scale_fill_manual(values = colores) +
    geom_jitter(
      aes(size = n_ind),
      alpha = 0.4,
      shape = 21,
      color = 'black',
      width = 0.4,
      height = 5) + 
    scale_y_continuous(breaks = seq(10, 140, by = 10),
                       limits = c(10, 140)) +
    theme_bw() +
    guides(fill = "none",
           color = "none",
           alpha = 'none') +
    labs(
      y = "AKDE 50% (ha)",
      title = "Fruit ripening",
      fill = "Study site",
      size = "Group Size") +
    theme(
      text = element_text(
        family = "Calibri",
        size = 34),
      axis.text.x = element_text(
        size = 0,
        angle = 55,
        color = "grey18",
        vjust = 0.6,
        hjust = 0.5),
      axis.title.x = element_text(size = 0),
      axis.text.y = element_text(color = 'grey15'),
      axis.title.y = element_text(vjust = 3),
      legend.box = element_blank(),
      plot.margin = unit(c(1, 1, 1, 1), "cm"))
*Figure 2: Predicted monthly Auto-correlated Kernel Density Estimator (AKDE) 50% (Core areas) of black capuchins (Sapajus nigritus) in pine plantations (orange) and Iguazú National Park (INP, green) across low and high ripe fruit availability seasons. Core areas were similar in both study sites across fruit ripening seasons*

Figure 2: Predicted monthly Auto-correlated Kernel Density Estimator (AKDE) 50% (Core areas) of black capuchins (Sapajus nigritus) in pine plantations (orange) and Iguazú National Park (INP, green) across low and high ripe fruit availability seasons. Core areas were similar in both study sites across fruit ripening seasons 

r_squared <- r.squaredGLMM(m3_50)
r2c <- r_squared[1, 2] 

set.seed(123)
folds <- createFolds(hr_month$log_akde_50_area, k = 10, list = TRUE)
errors <- numeric(length(folds))

for (i in seq_along(folds)) {
  train_data <- hr_month[-folds[[i]], ]  
  test_data <- hr_month[folds[[i]], ]   
  
  model <- lmer(log_akde_50_area ~ area + productividad + n_ind + (1 | grupo), 
                data = train_data)
  
  predictions <- predict(model, newdata = test_data, allow.new.levels = TRUE)
  
  errors[i] <- mean((predictions - test_data$log_akde_50_area)^2)
}

mean_cv_error <- mean(errors)
rmse_log <- sqrt(mean_cv_error)  
typical_error_percentage <- (exp(rmse_log) - 1) * 100  

model_performance <- data.frame(
  Metric = c("Marginal R² (R²m)", "Conditional R² (R²c)", "RMSE (log scale)", "Typical Error (%)"),
  Value = c(r2m, r2c, rmse_log, typical_error_percentage))

model_performance %>%
  flextable() %>%
  theme_tron_legacy() %>%  
  colformat_num(j = "Value", digits = 3) %>%
  bold(part = "header") %>%
  autofit() %>%
  set_caption("Model Performance Summary: R² and RMSE")
Model Performance Summary: R² and RMSE

Metric

Value

Marginal R² (R²m)

0.5722350

Conditional R² (R²c)

0.7725882

RMSE (log scale)

0.2913537

Typical Error (%)

33.8237844

5.3 Plantation Group - Bark stripping

hr_month_arauco2 <- readRDS ("files/hr_month_arauco3.rds")
akde50_bark_paper <-
  ggplot(hr_month_arauco2, aes(x = media_dano, y = akde_50_area)) +
    geom_jitter(
      size = 7.4, 
      alpha = 0.2, 
      color = 'black', 
      shape = 21, 
      fill = '#968282',
      width = 0.01,
      height = 5) +
    geom_smooth(
      method = "lm",
      linewidth = 3.5, 
      color = "#250303",
      fill = "#968282",
      se = TRUE,
      alpha = 0.2
    ) +
    labs(
      y = "Monthly AKDE 50% (ha)",
      x = expression("Monthly bark-stripping (m"^2*")")) +
    scale_y_continuous(breaks = seq(0, 140, by = 20), limits = c(0,140)) +
    scale_x_continuous(breaks = seq(0, 1.5, by = 0.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.1, vjust = -0.005),
      axis.text.y = element_text( color = "grey15"),
      axis.title.x = element_text( margin = margin(t = 10)),
      axis.title.y = element_text( vjust = 3),
      plot.margin = unit(c(1, 1, 1, 1), "cm"))


akde95_bark_paper <- 
  ggplot(hr_month_arauco2, aes(x = media_dano, y = akde_95_area)) +
  geom_jitter(
    size = 7.4, 
    alpha = 0.2, 
    color = 'black', 
    shape = 21, 
    fill = '#968282',
    width = 0.01,
    height = 5) +
  geom_smooth(
    method = "lm",
    linewidth = 3.5, 
    color = "#250303",
    fill = "#968282",
    se = TRUE,
    alpha = 0.2
  ) +
  labs(
    y = "Monthly AKDE 95% (ha)",
    x = expression("Monthly bark-stripping (m"^2*")")) +
#  scale_y_continuous(breaks = seq(0, 140, by = 20), limits = c(0,140)) +
  scale_x_continuous(breaks = seq(0, 1.5, by = 0.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.1, vjust = -0.005),
    axis.text.y = element_text(color = "grey15"),
    axis.title.x = element_text( margin = margin(t = 10)),
    axis.title.y = element_text( vjust = 3),
    plot.margin = unit(c(1, 1, 1, 1), "cm"))
  
akde95_bark_paper + akde50_bark_paper
*Figure S2. Monthly Auto-correlated Kernel Density Estimator (AKDE) 95% (home ranges; B) and 50% (core areas; C) of black capuchins (Sapajus nigritus) in pine plantations as a function of the average monthly bark-stripped area of pine trees.*

Figure S2. Monthly Auto-correlated Kernel Density Estimator (AKDE) 95% (home ranges; B) and 50% (core areas; C) of black capuchins (Sapajus nigritus) in pine plantations as a function of the average monthly bark-stripped area of pine trees.