Data and code for: Trade-offs across life history stages and social association types shape winter communal roosting in a long-lived raptor

Trade-offs across life history stages and social association types shape winter communal roosting in a long-lived raptor Code to generate communal roosting estimates using GPS-tracking data and model communal roosting propensities Table of contents 1. Identify individuals communal roosting . . . . . . . . . . . . . . . . . . . . . 3 2.Dataannotation.................................. 8 3.Summary...................................... 15 4.Dataformatting.................................. 15 5. Explore the relationship between roosting propensity across the season (in juvsandadults) ............................... 19 6. Model (roosting ~ age, sex, breeding status) . . . . . . . . . . . . . . . . . . 20 Juveniles........................................ 20 Adults......................................... 25 7. Disentangling elective disappearance vs behavioural plasticity . . . . . . . . 27 1 Figure 1: A Red kite roost in the Canton of Fribourg, Switzerland. Photo: Patrick Scherler # Clear R's mind rm(list = ls()) # Load libraries library(here) # easy referencing library(lubridate) # data wrangling library(dplyr) # data wrangling library(tidyverse) # data wrangling library(ggplot2) # plot library(sf) # spatial data wrangling library(data.table) # data wrangling library(brms) # model library(performance) # for the bayesian R2 library(bayestestR) # for the describe_posterior function library(tidybayes) # for the add_epred_draws function library(modelr) # to generate simulated dataset library(vcd) # for Cramers' coefficient 2 library(zoo) # data wrangling library(broom.mixed) # to export tidy model outputs library(readxl) 1. Identify individuals communal roosting We use nightly-averaged, individual GPS location to perform a spatial join to identify locations that intersect with each other. This will give you pairs of individuals that are within 300 meters of each other. The generated dataframe (“spatial_join”) consists of the pairs of individuals within 300m of each other. This is repeated for each temporal interval roost_dayweek <- list() for (j in 1:14) { results <- list() unique_nights <- na.omit(unique(nightmed[[paste0("interval",j)]])) for (i in 1:length(unique_nights)) { # Filter data for the current night data_night <- nightmed[which(nightmed[[paste0("interval",j)]] == unique_nights[i]),]↪ # Create a spatial buffer for the current night buffer_radius <- 300 # in meters data_night_buffer <- st_buffer(data_night, dist = buffer_radius) # Perform a spatial join for the current night spatial_join <- st_join(data_night_buffer, data_night, join = st_intersects)↪ # Exclude pairs where the same individual is present in both columns proximity_pairs <- subset(spatial_join, id.x != id.y) proximity_pairs$aggr <- as.character(proximity_pairs[[paste0("interval",j,".x")]])↪ # Add to the dataset also the individuals that don't intersect with anyone toadd <- data_night %>% filter(!id %in% proximity_pairs$id.x) 3 colnames(toadd) <- paste0(colnames(toadd), ".x") colnames(toadd)[which(colnames(toadd) == "geometry.x")] <- "geometry" toadd2 <- toadd colnames(toadd2) <- str_replace(colnames(toadd2),".x",".y") toadd2 <- toadd2 %>% st_drop_geometry() toadd <- cbind(toadd, toadd2) toadd$aggr <- as.character(toadd[[paste0("interval",j,".x")]]) toadd <- toadd[,c(order(colnames(proximity_pairs)))] toadd$id.y <- NA proximity_pairs <- rbind(proximity_pairs, toadd) # Store the results for the current night results[[i]] <- proximity_pairs } roost_dayweek[[j]] <- results } Associate the brood ID to each individual and year dt <- as.data.table(combined_results) id_brood <- unique(dt, by = c("id.x","winter.x")) %>% dplyr::rename(id = id.x) # combination of each unique individual per winter↪ # For individuals tagged as juveniles, until they start breeding, the brood_ID is the one of origin↪ attributes <- read_csv(here::here("./Data/brood_info.csv")) %>% dplyr::select(id, age, hatch_year, sex_compiled, nest_of_origin_ID, brood_ID, nest_ID_15, nest_ID_16, nest_ID_17, nest_ID_18, nest_ID_19, nest_ID_20, nest_ID_21, territoryID_2016, ↪ ↪ territoryID_2017, territoryID_2018, territoryID_2019, territoryID_2020, territoryID_2021 ,tag_year, manufacturer) # ids tagged as adults adults_attr <- attributes %>% 4 filter(!age == "1CY")%>% dplyr::select(id, nest_ID_15, nest_ID_16, nest_ID_17, nest_ID_18, nest_ID_19, nest_ID_20, nest_ID_21) %>%↪ mutate_all(as.character) %>% pivot_longer(!id, names_to = "winter.x",values_to = "nest_ID")%>%↪ mutate(winter.x = sub(".*_(\\d+)","20\\1", winter.x)) %>% # transform the nest year to brood year↪ mutate(brood_ID = paste(nest_ID, winter.x, sep=" "), age = "adult")↪ # ids tagged as juveniles - their brood of origin juvs_attr <- attributes %>% filter(age == "1CY")%>% dplyr::select(id, hatch_year, nest_of_origin_ID, brood_ID) %>% dplyr::rename(winter.x = hatch_year, nest_ID = nest_of_origin_ID)↪ # ids tagged as juveniles - their brood when they start breeding settl_attr <- attributes %>% filter(age == "1CY")%>% dplyr::select(id, nest_ID_15, nest_ID_16, nest_ID_17, nest_ID_18, nest_ID_19, nest_ID_20, nest_ID_21) %>%↪ mutate_all(as.character) %>% pivot_longer(!id, names_to = "winter.x",values_to = "nest_ID") %>%↪ mutate(winter.x = sub(".*_(\\d+)","20\\1", winter.x)) %>% # transform the nest year to brood year↪ mutate(brood_ID = paste(nest_ID, winter.x, sep=" ")) %>% drop_na(nest_ID) juvs <- rbind(juvs_attr, settl_attr) # for each juvenile, we add the complete history year-history. # The brood ID will first be the brood of origin and it switches when it starts breeding↪ juvs$winter.x <- as.numeric(juvs$winter.x) juvs <- juvs %>% arrange(id, winter.x) %>% group_by(id) %>% 5 complete(winter.x = full_seq(min(winter.x):2021,1)) %>% mutate(nest_ID = na.locf(nest_ID), brood_ID = na.locf(brood_ID)) %>% ungroup() %>% mutate(age = "juv") #put together all info all_attr <- rbind(adults_attr, juvs) # the NAs are the adults that stopped transmitting or died↪ all_attr$id_year <- paste(all_attr$id, all_attr$winter.x, sep="_")# to join it to the original observed dataset↪ # create a id_winter column and then join it to the main dataset (combined_results)↪ combined_results <- combined_results %>% mutate(id_year = paste(id.x, winter.x, sep="_"))↪ combined_results <- combined_results %>% left_join(all_attr[,c("brood_ID","id_year")] %>% dplyr::rename(brood_ID.x =brood_ID)) ↪ ↪ # add the brood ID of the second individual of the couple combined_results <- combined_results %>% mutate(id_year = paste(id.y, winter.x, sep="_"))↪ all_attr <- all_attr %>% dplyr::rename(brood_ID.y = brood_ID) # For the 3 individuals that are tagged as juveniles and then start breeding # we have to change the age (6, 322 and 91) # 6 from 2019 on # 322 and 91 from 2020 on all_attr$age[which(all_attr$id == "6" &all_attr$winter.x %in% c("2019", "2020","2021"))] <- "adult"↪ all_attr$age[which(all_attr$id == "322" &all_attr$winter.x %in% c("2020", "2021"))] <- "adult"↪ all_attr$age[which(all_attr$id == "91" &all_attr$winter.x %in% c("2020", "2021"))] <- "adult"↪ combined_results <- combined_results %>% left_join(all_attr[,c("brood_ID.y","id_year")], by = "id_year")↪ # if there is "unknown" in the brood ID is because the bird was caught without finding the nest --> we make it NA↪ combined_results <- combined_results %>% 6 mutate(brood_ID.x = ifelse(grepl("unknown", brood_ID.x, ignore.case = TRUE), NA, brood_ID.x),↪ brood_ID.y = ifelse(grepl("unknown", brood_ID.y, ignore.case = TRUE), NA, brood_ID.y))↪ combined_results$id_unique <- seq(1:length(combined_results$id.x)) # Export association data for the # parent-offspring and kin-kin association # analysis later # combined_results %>% # filter(temp_thresh == 5) %>% # write_rds(here::here("./Data/Association_types/rk_300_5night.rds")) # As well as the individual attributes # all_attr %>% # write_rds(here::here("./Data/Association_types/all_attr.rds")) Define an individual as roosting or not depending on the number of individuals it’s roosting with # How many other individuals were in the same location that night? If more than 2, define as communal roosting↪ dt <- as.data.table(combined_results) roosts_tot <- dt[, .(n_ind = n_distinct(id.y, na.rm=T)), by = .(temp_thresh, id.x, aggr)] %>%↪ as_tibble() %>% mutate(roosting = if_else(n_ind >= 2,1,0)) # attach the coordinates of the night locations to plot roosts_xy <- unique(dt, by = c("id.x","temp_thresh","aggr","doy.x")) roosts_tot <- roosts_tot %>% left_join(roosts_xy[,c("id.x","winter.x","temp_thresh", "aggr","doy.x","geometry","brood_ID.x")], by = c("id.x","temp_thresh","aggr")) %>% ↪ ↪ distinct(temp_thresh, aggr, id.x, .keep_all = T) %>% st_as_sf() rm(dt, roosts_xy) 7 # This roosts_tot file will be needed for the validation with the count data # (Supplementary Information) saveRDS(roosts_tot, here::here("./Data/roosts_tot.rds")) Based on sensitivity analysis and comparison with communal roosting census data (see Supp. Information), we decided to use 5 nights as temporal threshold and 300 m as spatial threshold 2. Data annotation Add individual information (age, breeding status, sex) roosts_tot_5 <- roosts_tot_5 %>% dplyr::rename(id = id.x) %>% dplyr::select(-c(temp_thresh)) %>% left_join(attributes[,c("id","hatch_year","sex_compiled","tag_year")], by = "id") ↪ ↪ # add age roosts_tot_5 <- roosts_tot_5 %>% group_by(id) %>% mutate(age = (as.numeric(winter.x) - as.numeric(hatch_year))+1)%>% # the NAs are adults. we add 1 because 0 is the 1st CY ↪ ↪ ungroup() %>% filter(!age %in% c(8:18)) roosts_tot_5$age[which(is.na(roosts_tot_5$age))] <- "3+" # add breeding status of the year before roosts_tot_5$id_year <- paste(roosts_tot_5$id, roosts_tot_5$winter.x, sep ="_")↪ roosts_4 <- list() for (i in 1:length(unique(roosts_tot_5$id_year))) { idx <- roosts_tot_5 %>% filter(id_year == unique(id_year)[i]) idx_territory <- attributes %>% filter(id == unique(idx$id)) %>% dplyr::select(paste0("territoryID_",unique(idx$winter.x))) idx$breeding <- as.numeric(idx_territory) 8 idx <- idx %>% mutate(breeding = if_else(is.na(breeding),"0","1")) roosts_4[[i]] <- idx } roosts_tot_5 <- do.call(bind_rows, roosts_4) rm(roosts_4) unique(roosts_tot_5$age) Add distance to the study area as weight # Add distance to the study area as a covariate lon_extent <- c(7.144554,7.373007)# longitudinal vertices of the study area lat_extent <- c(46.68084,46.90619)# latitudinal vertices of the study area # add distance to the margin instead of lat-lon bbox <- data.frame(lon = lon_extent, lat = lat_extent) %>% st_as_sf(coords = c("lon","lat"), crs = 4326)%>% st_transform(2056)%>% st_bbox() # Create a matrix representing the bounding box as a polygon bbox_matrix <- matrix(c(bbox['xmin'], bbox['ymin'], bbox['xmax'], bbox['ymin'], bbox['xmax'], bbox['ymax'], bbox['xmin'], bbox['ymax'], bbox['xmin'], bbox['ymin']), ncol = 2,byrow = TRUE) # Create an sf polygon object bbox <- st_polygon(list(bbox_matrix)) %>% st_sfc(crs = 2056) roosts_tot_5$sa <- as.numeric(st_within(roosts_tot_5, bbox)) # now it's either 0 (in the study area) or NA (not in the study area)↪ roosts_tot_5$sa[which(is.na(roosts_tot_5$sa))] <- "0" #visual check ggplot() + geom_sf(data = bbox, alpha = .2)+ geom_sf(data = roosts_tot_5, mapping = aes(col = sa)) + theme_bw() 9 roosting_juv$age.z <- scale(roosting_juv$age)[,1] roosting_juv <- roosting_juv %>% arrange(id,winter.x,doy.x) # adults roosting_ad$id <- as.factor(roosting_ad$id) roosting_ad$doy.x[which(roosting_ad$doy.x <100)] <- roosting_ad$doy.x[which(roosting_ad$doy.x <100)] +365↪ roosting_ad$doy.z <- scale(roosting_ad$doy.x, center = TRUE)[,1] roosting_ad$winter.x <- as.factor(roosting_ad$winter.x) roosting_ad$roosting <- as.factor(roosting_ad$roosting) roosting_ad$sex_compiled <- as.factor(roosting_ad$sex_compiled) roosting_ad <- roosting_ad %>% arrange(id,winter.x,doy.x) # What are the sample sizes per categories? roosting_juv %>% group_by(age, breeding, sex_compiled) %>% dplyr::summarise(n_obs = n(), n_ind = n_distinct(id)) # A tibble: 22 x 5 # Groups: age, breeding [11] age breeding sex_compiled n_obs n_ind <dbl> <fct> <fct> <int> <int> 1 1 0 f 121 14 2 1 0 m 173 19 3 2 0 f 149 17 4 2 0 m 491 45 5 3 0 f 201 19 6 3 0 m 545 53 7 3 1 f 141 13 8 3 1 m 45 4 9 4 0 f 100 8 10 4 0 m 279 20 # i 12 more rows roosting_ad %>% group_by(sex_compiled) %>% dplyr::summarise(n_obs = n(), n_ind = n_distinct(id)) # A tibble: 2 x 3 16 sex_compiled n_obs n_ind <fct> <int> <int> 1 f 1767 31 2 m 2010 37 17 18 5. Explore the relationship between roosting propensity across the season (in juvs and adults) 320 340 360 380 0.0 0.2 0.4 0.6 0.8 JUVENILES doy proportion of nights spent at communal roosts 320 340 360 380 0.0 0.2 0.4 0.6 0.8 1.0 ADULTS doy proportion of nights spent at communal roosts 19 There is a non-significant negative relationship between day of the year and the probability or roosting communally in juveniles and positive significant relationship in individuals tagged as adults: Person’s r for juveniles: -0.148, p-value: 0.16 Person’s r for adults: 0.266, p-value: 0.01 6. Model (roosting ~ age, sex, breeding status) Juveniles # Create the numerical day sequence roosting_juv <- roosting_juv %>% mutate(aggr = as.numeric(sub(".*_","", roosting_juv$aggr))) %>%↪ group_by(id, winter.x) %>% mutate(doy_num = aggr -min(aggr) +1)%>% ungroup() %>% mutate(doy_num.z = scale(doy_num)[,1]) # Add weights that go from 0 to 1 roosting_juv$weights <- scale(roosting_juv$sa, center = min(roosting_juv$sa), scale = max(roosting_juv$sa) -min(roosting_juv$sa))↪ roosting_juv$weights <- as.numeric(1-roosting_juv$weights) range(roosting_juv$weights) [1] 0 1 roosting_full_quadr <- brm(bf(roosting|weights(weights) ~age.z +I(age.z^2) +sex_compiled +↪ breeding +(doy_num.z|id_year) +(1|id) + (1|winter.x)),↪ data=roosting_juv, family=bernoulli(), iter = 4000, cores = 4) # Export the table of the final model describe_posterior(roosting_full_quadr, effects="all")%>% mutate_if(is.numeric, function(x) round(x, digits = 3)) %>% dplyr::select(1:6)%>% 20 write_csv("./Results/Tables/Mod_juv.csv") ĎVisualisation Model check # 1. Binned residual plots: The following shows three binned residual plots, with each point showing y1−y, where y1 is based on simulated data from the posterior predictive distributions, and y is the observed data. Note that we need the binned residuals or predictive errors, as the prediction error is either 0 or 1, as shown in the Figure below. the binned margins were based on the observed data, whereas the dots were predictive errors from replicated data. Predictive proportions (x-axis) ↪ ↪ ↪ ↪ ↪ ↪ pp_check(roosting_full_quadr, type = "error_binned",nsamples = 9) 0.00 0.25 0.50 0.75 1.00 0.00 0.25 0.50 0.75 1.00 0.00 0.25 0.50 0.75 1.00 −0.2 −0.1 0.0 0.1 0.2 −0.2 −0.1 0.0 0.1 0.2 −0.2 −0.1 0.0 0.1 0.2 Predicted proportion Average Errors (with 2SE bounds) # 2. Average classification accuracy --> OBSERVED VS SIMULATED: Another way to evaluate a logistic regression model is to look at classification error, which is analogous to R^2 for normal regression.The simplest measure is to assign observations with predicted probabilities larger than �0 to have a value of 1, and to assign observations with predicted probabilities smaller than �0 to have a value of 0, where �0 is some chosen cutoff, and is usually chosen as the proportion of 1s in the sample. For example, for our model, ↪ ↪ ↪ ↪ ↪ ↪ ↪ 21 m1_pred <- predict(roosting_full_quadr, type = "response")[ , "Estimate"] m1_pred <- as.numeric(m1_pred > mean(as.numeric(as.character(roosting_juv$roosting))))↪ # Classification table (classtab_m1 <- table(predicted = m1_pred, observed = roosting_juv$roosting)) observed predicted 0 1 0 1863 231 1 312 1746 # Average classification accuracy (acc_m1 <- sum(diag(classtab_m1)) /sum(classtab_m1)) [1] 0.8692197 Plot the age and sex effect Age_sex_effect <- roosting_juv %>% data_grid(sex_compiled, age.z, breeding) %>% add_epred_draws(roosting_full_quadr, re_formula=NA,ndraws =100)%>%↪ mutate(age = age.z*sd(roosting_juv$age)+mean(roosting_juv$age))↪ Age_sex_effect_b <- Age_sex_effect[which(Age_sex_effect$breeding == "1"),] Age_sex_effect_nb <- Age_sex_effect[which(Age_sex_effect$breeding == "0"),] # rawdata for the points rawdata <- roosting_juv rawdata$roosting <- as.numeric(as.character(rawdata$roosting)) rawdata$roosting[which(rawdata$sex_compiled == "f"&rawdata$roosting == 0)] <- -0.05↪ rawdata$roosting[which(rawdata$sex_compiled == "m"&rawdata$roosting == 1)] <- 1.05↪ # Plot only the non-breeders effect_grouped_nb <- Age_sex_effect_nb %>% 22 group_by(age, sex_compiled) %>% dplyr::summarise(med_pred = median(.epred)) %>% ungroup effect_grouped <- Age_sex_effect %>% group_by(age, sex_compiled, breeding) %>% dplyr::summarise(med_pred = median(.epred)) %>% ungroup # Remove from the figure the breeders at age 1 (there are none in reality) Age_sex_effect$.epred[which(Age_sex_effect$breeding == 1&Age_sex_effect$age == 1)] <- NA↪ effect_grouped$med_pred[which(effect_grouped$breeding == 1& effect_grouped$age == 1)] <- NA↪ Age_sex_effect$.epred[which(Age_sex_effect$breeding == 1&Age_sex_effect$age == 2)] <- NA↪ effect_grouped$med_pred[which(effect_grouped$breeding == 1& effect_grouped$age == 2)] <- NA↪ fig2 <- ggplot(Age_sex_effect, aes(y = .epred, x = age, fill=sex_compiled, group = interaction(sex_compiled, breeding))) + stat_lineribbon(alpha=.5,.width=.95)+ geom_point(data = effect_grouped, aes(y = med_pred,x = age, group = interaction(sex_compiled, breeding)),↪ size = 3.5)+ scale_color_manual(values=c("khaki3","steelblue4")) + scale_fill_manual(values=c("khaki3","steelblue4")) + xlab("Age [years]")+ylab("% Time spent in a communal roost")+↪ scale_x_continuous(breaks=seq(1,7,1))+ scale_y_continuous(breaks=seq(0,1,0.2)) + facet_grid(.~breeding) + theme_classic()+ theme(legend.position = "none", axis.text = element_text(size = 15), axis.title = element_text(size = 15), strip.background = element_blank(), # remove facet_grid labels↪ strip.text.x = element_blank()) 23 fig2 12345671234567 0.0 0.2 0.4 0.6 0.8 Age [years] % Time spent in a communal roost # ggsave(fig2, file = "./Results/Plots/Fig2.png", # width = 8, height = 5) # Extract numbers for the manuscript - median and 95 % CrI for # Age categories Age_sex_effect %>% filter(sex_compiled == levels(sex_compiled)[2], breeding == levels(breeding)[1]) %>%↪ group_by(age) %>% dplyr::summarise(m_prob = median(.epred), low_cri = quantile(.epred, 0.025), up_cri = quantile(.epred, 0.975)) %>% mutate(across(c(m_prob, low_cri, up_cri), ~round(., 2))) # A tibble: 7 x 4 age m_prob low_cri up_cri <dbl> <dbl> <dbl> <dbl> 1 1 0.56 0.28 0.81 2 2 0.74 0.6 0.86 24 3 3 0.79 0.65 0.89 4 4 0.76 0.58 0.88 5 5 0.64 0.42 0.83 6 6 0.37 0.15 0.73 7 7 0.1 0.02 0.48 # Sex Age_sex_effect %>% filter(breeding == levels(breeding)[1], age == 4)%>% group_by(sex_compiled) %>% summarise(m_prob = median(.epred), low_cri = quantile(.epred, 0.025), up_cri = quantile(.epred, 0.975)) %>% mutate(across(c(m_prob, low_cri, up_cri), ~round(., 2))) # A tibble: 2 x 4 sex_compiled m_prob low_cri up_cri <fct> <dbl> <dbl> <dbl> 1 f 0.37 0.18 0.59 2 m 0.76 0.58 0.88 # Breeding effect Age_sex_effect %>% filter(age == 4, sex_compiled == levels(sex_compiled)[1]) %>% group_by(breeding) %>% summarise(m_prob = median(.epred), low_cri = quantile(.epred, 0.025), up_cri = quantile(.epred, 0.975)) %>% mutate(across(c(m_prob, low_cri, up_cri), ~round(., 2))) # A tibble: 2 x 4 breeding m_prob low_cri up_cri <fct> <dbl> <dbl> <dbl> 1 0 0.37 0.18 0.59 2 1 0.14 0.07 0.27 Adults 25