Test different interaction methods
Johannes Krietsch
2023-08-15
#==============================================================================================================Data and code from “Mutual mate guarding and limited sexual conflict
in a sex-role reversed shorebird” Contributor: Johannes Krietsch
📍 This script runs relative to the project’s root directory and
describes how I tested different methods to find the best way to define
if birds were “together”.
#==============================================================================================================
# Packages
sapply( c('data.table', 'magrittr', 'sdb', 'ggplot2', 'anytime', 'viridis', 'auksRuak', 'foreach', 'sf', 'knitr',
'stringr', 'windR', 'ggnewscale', 'doFuture', 'patchwork'),
require, character.only = TRUE)## data.table magrittr sdb ggplot2 anytime viridis auksRuak foreach sf knitr stringr windR ggnewscale
## TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## doFuture patchwork
## TRUE TRUE# Functions
source('./R/0_functions.R')
# Lines to run to create html output
opts_knit$set(root.dir = rprojroot::find_rstudio_root_file())
# rmarkdown::render('./R/4_interaction_method.R', output_dir = './OUTPUTS/R_COMPILED')
# Projection
PROJ = '+proj=laea +lat_0=90 +lon_0=-156.653428 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 '
# Data
d = fread('./DATA/NANO_TAGS.txt', sep = '\t', header = TRUE, nThread = 20) %>% data.table
d = d[filtered == TRUE]
st_transform_DT(d)
dp = fread('./DATA/PRODUCTS/PAIR_WISE_DIST_CLOSEST.txt', sep = '\t', header = TRUE, nThread = 20) %>% data.table
dn = fread('./DATA/NESTS.txt', sep = '\t', header = TRUE, nThread = 20) %>% data.table
dn[, initiation_y := as.POSIXct(format(initiation, format = '%m-%d %H:%M:%S'), format = '%m-%d %H:%M:%S', tz = 'UTC')]
# change projection
st_transform_DT(dn)Define breeding pairs with both sexes tagged
# start and end of the data
d[, start := min(datetime_), by = ID]
d[, end := max(datetime_), by = ID]
dID = unique(d, by = 'ID')
# check if data overlap
dn = merge(dn, dID[, .(male_id = ID, start_m = start, end_m = end)], by = 'male_id', all.x = TRUE)
dn = merge(dn, dID[, .(female_id = ID, start_f = start, end_f = end)], by = 'female_id', all.x = TRUE)
# subset nests with both IDs tagged
dn = dn[!is.na(start_m) & !is.na(start_f)]
# subset nests with both IDs tagged and overlapping time intervals
dn[, overlap := DescTools::Overlap(c(start_m, end_m), c(start_f, end_f)), by = nestID]
dn = dn[overlap > 0]
# check overlap with initiation date
dn[, overlap_initiation_m := DescTools::Overlap(c(start_m, end_m), c(initiation - 86400, initiation + 86400)),
by = nestID]
dn[, overlap_initiation_f := DescTools::Overlap(c(start_f, end_f), c(initiation - 86400, initiation + 86400)),
by = nestID]
dn = dn[overlap_initiation_m > 0 & overlap_initiation_f > 0]
# nest data
dnID = dn[, .(year_, nestID, male_id, female_id, initiation, initiation_y, lat_n = lat, lon_n = lon)]
dnID = unique(dnID, by = 'nestID')
# as integer
dnID[, male_id := as.integer(male_id)]
dnID[, female_id := as.integer(female_id)]
# # create directory for each of these breeding pairs
# dnID[, directory := paste0('//ds/grpkempenaers/Hannes/temp/PAIRS_ANIMATION_EACH/', nestID)]
# dnID[, dir.create(file.path(directory), showWarnings = FALSE), by = 1:nrow(dnID)]
# unique IDs
IDu = unique(c(dnID[, male_id], dnID[, female_id]))
# subset d
d = d[ID %in% IDu]
# merge with nest and initiation date
dnIDu = rbind(dnID[, .(year_, ID = female_id, nestID, initiation, sex = 'F')],
dnID[, .(year_, ID = male_id, nestID, initiation, sex = 'M')])
d = merge(d, dnIDu[, .(year_, ID, nestID, initiation, sex)], by = c('ID', 'year_'), all.x = TRUE,
allow.cartesian = TRUE)
d = d[!is.na(nestID)]Define interactions
# merge with nests
dp = merge(dp, dnID, by.x = c('ID1', 'ID2'), by.y = c('male_id', 'female_id'))
# interactions
dp[, interaction := distance_pair < 30]
# count bouts of split and merge
dp[, bout := bCounter(interaction), by = nestID]
dp[, bout_seq := seq_len(.N), by = .(nestID, bout)]
dp[, bout_seq_max := max(bout_seq), by = .(nestID, bout)]
dp[interaction == FALSE & bout_seq_max == 1, interaction := TRUE]
# interaction ID
dp[, int_id := seq_len(.N), by = nestID]
# first and last interaction
dp[interaction == TRUE, first_int := min(datetime_1), by = nestID]
dp[, first_int := min(first_int, na.rm = TRUE), by = nestID]
dp[interaction == TRUE, last_int := max(datetime_1), by = nestID]
dp[, last_int := min(last_int, na.rm = TRUE), by = nestID]
# relative nest initiation date
dp[, initiation_rel := difftime(datetime_1, initiation, units = 'days') %>% as.numeric()]Plot for each nest
# sort by most time together before initiation
ds = dp[initiation_rel < 0 & interaction == TRUE, .N, by = nestID]
setorder(ds, -N)
# order nest ID
dp[, nestID := factor(nestID %>% as.factor, levels = ds[, nestID])]
ggplot(data = dp) +
geom_tile(aes(initiation_rel, nestID, fill = interaction), width = 0.5, show.legend = FALSE) +
scale_fill_manual(values = c('TRUE' = 'green4', 'FALSE' = 'firebrick3', 'NA' = 'grey50')) +
geom_vline(aes(xintercept = 0), color = 'black', size = 3, alpha = 0.5) +
geom_vline(aes(xintercept = -2), color = 'black', size = 3, alpha = 0.5) +
xlab('Date relative to initiation') + ylab('Nest') +
theme_classic()## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
Interactions using different distance thresholds
# distance threshold
threshold = sequence(15, 10, 10)
o = foreach(i = 1:length(threshold), .combine = 'rbind') %do% {
distance_threshold = threshold[i]
# interactions
dp[, interaction := distance_pair < distance_threshold]
# count bouts of split and merge
dp[, bout := bCounter(interaction), by = nestID]
dp[, bout_seq := seq_len(.N), by = .(nestID, bout)]
dp[, bout_seq_max := max(bout_seq), by = .(nestID, bout)]
dp[interaction == FALSE & bout_seq_max == 1, interaction := TRUE]
# round to days
dp[, initiation_rel := round(initiation_rel, 0)]
# daily points of both individuals
dp[, N_daily := .N, by = .(nestID, initiation_rel)]
# daily interactions
dp[interaction == TRUE, N_together := .N, by = .(nestID, initiation_rel)]
dp[interaction == FALSE, N_together := NA]
dp[, N_together := mean(N_together, na.rm = TRUE), by = .(nestID, initiation_rel)]
dp[is.na(N_together), N_together := 0]
# unique data
ds = unique(dp, by = c('nestID', 'initiation_rel'))
ds[, per_together := N_together / N_daily * 100]
# nests to exclude
n2 = c('R201_19', 'R231_19', 'R905_19', 'R502_19')
ds = ds[!(nestID %in% n2)]
# exclude pairs before mate guarding started
ds = ds[!(initiation_rel < 0 & per_together < 50)]
ds = ds[, .(per_together = median(per_together)), by = initiation_rel]
ds[, distance_threshold := distance_threshold]
ds
}
setorder(o, initiation_rel)
ggplot(data = o) +
geom_point(aes(initiation_rel, per_together, color = distance_threshold, group = distance_threshold),
size = 2, alpha = 1) +
geom_path(aes(initiation_rel, per_together, color = distance_threshold, group = distance_threshold),
size = 1, alpha = 0.5) +
scale_color_viridis(direction = -1, limits = c(10, 150), name = 'distance threshold') +
geom_vline(aes(xintercept = 0), color = 'firebrick2', size = 3, alpha = 0.3) +
geom_vline(aes(xintercept = 3), color = 'firebrick2', size = 1, alpha = 0.3) +
xlab('Day relative to clutch initiation (= 0)') + ylab('Percentage of positions together') +
theme_classic(base_size = 8) +
theme(legend.position = c(0.8, 0.8))
Interactions using different maximal bout length
# distance threshold
bout_seq_max_value = sequence(6, 0, 1)
o = foreach(i = 1:length(bout_seq_max_value), .combine = 'rbind') %do% {
bsm = bout_seq_max_value[i]
# interactions
dp[, interaction := distance_pair < 30]
# count bouts of split and merge
dp[, bout := bCounter(interaction), by = nestID]
dp[, bout_seq := seq_len(.N), by = .(nestID, bout)]
dp[, bout_seq_max := max(bout_seq), by = .(nestID, bout)]
dp[interaction == FALSE & bout_seq_max <= i, interaction := TRUE]
# round to days
dp[, initiation_rel := round(initiation_rel, 0)]
# daily points of both individuals
dp[, N_daily := .N, by = .(nestID, initiation_rel)]
# daily interactions
dp[interaction == TRUE, N_together := .N, by = .(nestID, initiation_rel)]
dp[interaction == FALSE, N_together := NA]
dp[, N_together := mean(N_together, na.rm = TRUE), by = .(nestID, initiation_rel)]
dp[is.na(N_together), N_together := 0]
# unique data
ds = unique(dp, by = c('nestID', 'initiation_rel'))
ds[, per_together := N_together / N_daily * 100]
# nests to exclude
n2 = c('R201_19', 'R231_19', 'R905_19', 'R502_19')
ds = ds[!(nestID %in% n2)]
# exclude pairs before mate guarding started
ds = ds[!(initiation_rel < 0 & per_together < 50)]
ds = ds[, .(per_together = median(per_together)), by = initiation_rel]
ds[, bout_seq_max_value := bsm]
ds
}
setorder(o, initiation_rel)
ggplot(data = o) +
geom_point(aes(initiation_rel, per_together, color = bout_seq_max_value, group = bout_seq_max_value),
size = 2, alpha = 1) +
geom_path(aes(initiation_rel, per_together, color = bout_seq_max_value, group = bout_seq_max_value),
size = 1, alpha = 0.5) +
scale_color_viridis(direction = -1, limits = c(0, 5), name = 'min bout lenght') +
geom_vline(aes(xintercept = 0), color = 'firebrick2', size = 3, alpha = 0.3) +
geom_vline(aes(xintercept = 3), color = 'firebrick2', size = 1, alpha = 0.3) +
xlab('Day relative to clutch initiation (= 0)') + ylab('Percentage of positions together') +
theme_classic(base_size = 8) +
theme(legend.position = c(0.8, 0.8))
Include movement before
dp = copy(dpc)
### positions before and after
# shift positions
dp[, lat1_before := shift(lat1, type = 'lag'), by = nestID]
dp[, lon1_before := shift(lon1, type = 'lag'), by = nestID]
dp[, lat2_before := shift(lat2, type = 'lag'), by = nestID]
dp[, lon2_before := shift(lon2, type = 'lag'), by = nestID]
dp[, lat1_next := shift(lat1, type = 'lead'), by = nestID]
dp[, lon1_next := shift(lon1, type = 'lead'), by = nestID]
dp[, lat2_next := shift(lat2, type = 'lead'), by = nestID]
dp[, lon2_next := shift(lon2, type = 'lead'), by = nestID]
# distance to position before and after
dp[, distance1_before := sqrt(sum((c(lon1, lat1) - c(lon1_before, lat1_before))^2)) , by = 1:nrow(dp)]
dp[, distance1_next := sqrt(sum((c(lon1, lat1) - c(lon1_next, lat1_next))^2)) , by = 1:nrow(dp)]
dp[, distance2_before := sqrt(sum((c(lon2, lat2) - c(lon2_before, lat2_before))^2)) , by = 1:nrow(dp)]
dp[, distance2_next := sqrt(sum((c(lon2, lat2) - c(lon2_next, lat2_next))^2)) , by = 1:nrow(dp)]
# interaction with "movement buffer"
dp[, interaction_buffer := distance_pair < c(max(distance1_before, distance2_before) + 30), by = 1:nrow(dp)]
# simple interactions
dp[, interaction_threshold := distance_pair < distance_threshold]
dp[is.na(interaction), interaction := interaction_threshold]
# count bouts of split and merge
dp[, bout := bCounter(interaction), by = nestID]
dp[, bout_seq := seq_len(.N), by = .(nestID, bout)]
dp[, bout_seq_max := max(bout_seq), by = .(nestID, bout)]
dp[, any_interaction_threshold := any(interaction_threshold == TRUE), by = .(nestID, bout)]
dp[any_interaction_threshold == FALSE, interaction := FALSE]
# split points and merging points
dp[, interaction_next := shift(interaction, type = 'lead'), by = nestID]
dp[, interaction_before := shift(interaction, type = 'lag'), by = nestID]
# correct for true splits
dp[interaction == TRUE & interaction_next == FALSE & distance_pair > distance_threshold, interaction := FALSE]
# datetime relative to nest initiation date
dp[, datetime_rel_pair := difftime(datetime_1, initiation, units = 'days') %>% as.numeric()]
dp[, datetime_rel_pair0 := round(datetime_rel_pair, 0)]
# with buffer
ggplot(data = dp) +
geom_tile(aes(initiation_rel, nestID, fill = interaction_buffer), width = 0.5, show.legend = FALSE) +
scale_fill_manual(values = c('TRUE' = 'green4', 'FALSE' = 'firebrick3', 'NA' = 'grey50')) +
geom_vline(aes(xintercept = 0), color = 'black', size = 3, alpha = 0.5) +
geom_vline(aes(xintercept = -2), color = 'black', size = 3, alpha = 0.5) +
xlab('Date relative to initiation') + ylab('Nest') +
theme_classic()
# without buffer
ggplot(data = dp) +
geom_tile(aes(initiation_rel, nestID, fill = interaction), width = 0.5, show.legend = FALSE) +
scale_fill_manual(values = c('TRUE' = 'green4', 'FALSE' = 'firebrick3', 'NA' = 'grey50')) +
geom_vline(aes(xintercept = 0), color = 'black', size = 3, alpha = 0.5) +
geom_vline(aes(xintercept = -2), color = 'black', size = 3, alpha = 0.5) +
xlab('Date relative to initiation') + ylab('Nest') +
theme_classic()
## R version 4.2.2 (2022-10-31 ucrt)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 19045)
##
## Matrix products: default
##
## locale:
## [1] LC_COLLATE=German_Germany.utf8 LC_CTYPE=German_Germany.utf8 LC_MONETARY=German_Germany.utf8 LC_NUMERIC=C
## [5] LC_TIME=German_Germany.utf8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] patchwork_1.1.2 doFuture_0.12.2 future_1.31.0 ggnewscale_0.4.8 windR_1.0.0 stringr_1.5.0 knitr_1.42 sf_1.0-13
## [9] foreach_1.5.2 auksRuak_0.1.0 viridis_0.6.2 viridisLite_0.4.1 anytime_0.3.9 ggplot2_3.4.0 sdb_2019.12 magrittr_2.0.3
## [17] data.table_1.14.8
##
## loaded via a namespace (and not attached):
## [1] sfext_0.1.0 lubridate_1.9.2 bit64_4.0.5 httr_1.4.4 rprojroot_2.0.3 bslib_0.4.2 tools_4.2.2
## [8] utf8_1.2.3 R6_2.5.1 KernSmooth_2.23-20 rgeos_0.6-1 DBI_1.1.3 colorspace_2.1-0 raster_3.6-14
## [15] withr_2.5.0 sp_1.6-0 tidyselect_1.2.0 gridExtra_2.3 Exact_3.2 RMySQL_0.10.25 bit_4.0.5
## [22] compiler_4.2.2 cli_3.6.0 expm_0.999-7 labeling_0.4.2 sass_0.4.5 scales_1.2.1 classInt_0.4-9
## [29] mvtnorm_1.1-3 proxy_0.4-27 askpass_1.1 digest_0.6.31 rmarkdown_2.20 pkgconfig_2.0.3 htmltools_0.5.4
## [36] parallelly_1.34.0 filenamr_0.1.0 highr_0.10 fastmap_1.1.0 rlang_1.0.6 RMariaDB_1.2.2 readxl_1.4.1
## [43] rstudioapi_0.14 farver_2.1.1 jquerylib_0.1.4 generics_0.1.3 jsonlite_1.8.4 dplyr_1.1.0 Matrix_1.5-1
## [50] Rcpp_1.0.10 DescTools_0.99.47 munsell_0.5.0 fansi_1.0.4 lifecycle_1.0.3 terra_1.7-3 stringi_1.7.12
## [57] yaml_2.3.7 snakecase_0.11.0 MASS_7.3-58.1 rootSolve_1.8.2.3 grid_4.2.2 parallel_4.2.2 listenv_0.9.0
## [64] lmom_2.9 lattice_0.20-45 hms_1.1.3 pillar_1.8.1 boot_1.3-28 gld_2.6.6 codetools_0.2-18
## [71] glue_1.6.2 evaluate_0.20 vctrs_0.5.2 cellranger_1.1.0 gtable_0.3.1 purrr_1.0.1 cachem_1.0.6
## [78] xfun_0.40 janitor_2.1.0 e1071_1.7-13 class_7.3-20 ggspatial_1.1.7 tibble_3.1.8 iterators_1.0.14
## [85] cliExtras_0.1.0 units_0.8-2 timechange_0.2.0 globals_0.16.2