# get number of individuals per groupn_per_group <-sapply(unique(group_traj$group), USE.NAMES =FALSE,function(x) {# get those for that group X <- group_traj[group_traj$group == x, ]# count number of individuals (No NAs) n_indiv <- (sum(!sapply(X, anyNA)) -3)/3return(n_indiv) })table(n_per_group)
n_per_group
2 3 4 5 6
7 8 7 8 4
1.2 Compute pairwise distances for natural and artificial groups
Code
pairwise_dists_list <-sapply(split(group_traj, group_traj$group), pairwise_dist, simplify =FALSE)pairwise_dists <-do.call(rbind, pairwise_dists_list)# remove those with unlikely distancepairwise_dists <- pairwise_dists[pairwise_dists$distance >0.1& pairwise_dists$distance <7, ]
1.3 Generate simulated group flight trajectories
Code
## split into single individual trajectories split in a listsplit_group_traj <-split(group_traj, group_traj$group)# get single individual trajectoresindiv_traj_list <-lapply(split_group_traj, function(x) { x <- x[, sapply(x, function(x) sum(is.na(x))) <10]# get number of individualssuppressWarnings(indv_num <-unique(na.omit(as.numeric(gsub(".*?([0-9]+).*","\\1", names(x)))))) indiv_trajs <-lapply(indv_num, function(y) x[, grep(y, names(x))])return(indiv_trajs)})indiv_traj_list <-unlist(indiv_traj_list, recursive =FALSE, use.names =TRUE)# remove last character of the list names so all have the# original group namenames(indiv_traj_list) <-gsub(".$", "", names(indiv_traj_list))# mean number of individuals per group size in original datamean_n_indiv <-floor(mean(table(n_per_group)))group_sizes <-as.numeric(names(table(n_per_group)))group_sizes <-rep(group_sizes, each = mean_n_indiv)sim_dists <-lapply(1:30, function(x) sim_group_traj(group_sizes,seed = x))sim_trajs_list <-lapply(1:30, function(x) sim_group_traj(group_sizes,seed = x, trajectories =TRUE))# remove those with unlikely distancesim_dists <-lapply(sim_dists, function(x) { x <- x[x$distance >0.1& x$distance <7.81, ]return(x)})
1.4 Describe data
Mean distance between individuals in a group flight:
Code
sim_dists_df <-do.call(rbind, sim_dists)# mean distance by group typeagg <-aggregate(distance ~ type, sim_dists_df, mean)agg$sd <-aggregate(distance ~ type, sim_dists_df, sd)$distance# mean distance by group typeagg2 <-aggregate(distance ~ type, pairwise_dists, mean)agg2$sd <-aggregate(distance ~ type, pairwise_dists, sd)$distancerbind(agg, agg2)
type
distance
sd
simulated
2.050082
0.5696759
artificial
2.005765
0.6803455
natural
1.670052
0.7256127
Closest distance between individuals in a group flight:
Code
sim_dists_df <-do.call(rbind, sim_dists)# mean distance by group typeagg_min <-aggregate(distance ~ type + group, sim_dists_df, min)agg <-aggregate(distance ~ type, agg_min, mean)agg$sd <-aggregate(distance ~ type, agg_min, sd)$distance# mean distance by group typeagg_min2 <-aggregate(distance ~ type + group, pairwise_dists, min)agg2 <-aggregate(distance ~ type, agg_min2, mean)agg2$sd <-aggregate(distance ~ type, agg_min2, sd)$distancerbind(agg, agg2)
type
distance
sd
simulated
0.6592624
0.3312211
artificial
0.9213170
0.5873041
natural
0.7190676
0.3595199
1.5 Plot raw data
1.5.1 Reduce dimensionality
1.6 Generate simulated group flight trajectories
Code
## split into single individual trajectories# split in a listsplit_group_traj <-split(group_traj, group_traj$group)# pca funfun1 <-function(x) prcomp(x, center =FALSE, scale. =FALSE)$x[,1]# horizontal plane fun# fun1 <- function(x) x[,3]indiv_traj_pc1 <-lapply(indiv_traj_list, fun1)pca_list <-lapply(unique(names(indiv_traj_pc1)), function(x){ pca_df <-data.frame(indiv_traj_pc1[names(indiv_traj_pc1) == x])names(pca_df) <-paste0("ind", 1:ncol(pca_df)) stacked_pca_df <-stack(pca_df) stacked_pca_df$frame <-1:nrow(pca_df) stacked_pca_df$group <- x stacked_pca_df$group.size <-ncol(pca_df)return(stacked_pca_df)})pca_df <-do.call(rbind, pca_list)# for simulated datasim_traj_pc1 <-lapply(sim_trajs_list[[1]], function(x) x[,3])# pca_sim_list <- lapply(unique, function(x){pca_sim_list <-lapply(paste("sim", 1:30, "", sep ="-"), function(x){ pca_df <-data.frame(sim_traj_pc1[grep(pattern = x, names(sim_traj_pc1))])names(pca_df) <-paste0("ind", 1:ncol(pca_df)) stacked_pca_df <-stack(pca_df) stacked_pca_df$frame <-1:nrow(pca_df) stacked_pca_df$group <- x stacked_pca_df$group.size <-ncol(pca_df)return(stacked_pca_df)})pca_sim_df <-do.call(rbind, pca_sim_list)pca_sim_df$type <-"simulated"pca_df$type <-sapply(pca_df$group, function(x) pairwise_dists$type[group_traj$group == x][1])#bindpca_df <-rbind(pca_df, pca_sim_df)# order groups by typepca_df$group <-factor(pca_df$group, levels =unique(pca_df$group[order(pca_df$type)]))pca_df <- pca_df[complete.cases(pca_df), ]pca_df$type <-factor(pca_df$type, levels =c("natural", "artificial", "simulated"))# table(pca_df$type, pca_df$group.size)#get 3 names of group for each type levelset.seed(1)selected_groups <-lapply(unique(pca_df$type), function(x) {unique(pca_df$group[pca_df$type == x & pca_df$group.size ==4])[1:3]})pca_df <- pca_df[pca_df$frame <26, ]ggplot(pca_df[pca_df$group %in%unlist(selected_groups), ], aes(x = frame, y = values, group = ind, color = type)) +# geom_point() + geom_smooth(se =FALSE, span =0.3) +theme_classic() +# theme(legend.position = "none") + labs(x ="Frame", y ="Horizontal plane\nposition (m)") +# ylim(c = c(-3, 3)) +facet_wrap( dir ="v", type~ group, #scales = "free_y",ncol =3) +theme(strip.background =element_blank(),strip.text.x =element_blank()) +scale_color_viridis_d(option ="G",end =0.8,alpha =0.6)
1.6.1 Raw distances
All group sizes combined:
Code
# add original data to simulated data and calculate proximitynat_sim_dist_list <-lapply(sim_dists, function(x) {return(rbind(x, pairwise_dists))})# remove distances above 10 mnat_sim_dist_list <-lapply(nat_sim_dist_list, function(x) { x <- x[x$distance <10, ] x$type <-factor(x$type, levels =c("natural", "artificial", "simulated"))return(x)})dat_all <- dat <- nat_sim_dist_list[[1]]dat_all$group.size <-"All sizes"dat <-rbind(dat, dat_all)dat$group.size.factor <-factor(dat$group.size, levels =c(1:10, "All sizes"))n_grps <-aggregate(group ~ group.size.factor + type, dat, function(x)length(unique(x)))n_grps$frames <-aggregate(group ~ group.size.factor + type, dat, function(x)length(x))$groupn_grps$indivs <- n_grps$group *as.numeric(as.character(n_grps$group.size.factor))n_grps$group.size.factor <-factor(n_grps$group.size.factor)n_grps$distance <-ifelse(n_grps$group.size.factor %in%2:3,-1.2,-0.9)n_grps$distance[1] <--0.8n_grps$indivs[n_grps$group.size.factor =="All sizes"] <-aggregate(indivs ~ type, n_grps, sum)$indivsn_grps$n.labels <-paste0(n_grps$group, " / ", n_grps$indivs, " / ", n_grps$frames)n_grps$n.labels[1] <-" n = 3 groups /\n 6 individuals /\n 149 frames"n_grps$n.labels[6] <-" n = 18 groups /\n 74 individuals /\n 888 frames"n_grps$distance[6] <--0.3n_grps$distance[c(11, 17)] <--0.4custom_labels <-unique(nat_sim_dist_list[[1]]$group.size)custom_labels <-paste(c("Group size = ", rep("", 5)), custom_labels)names(custom_labels) <-unique(nat_sim_dist_list[[1]]$group.size)custom_labels <-unique(dat$group.size)custom_labels <-paste(c("Group size = ", rep("", 5)), custom_labels)names(custom_labels) <-unique(dat$group.size)# mean distance by group typeagg_gp <-aggregate(distance ~ type + group + group.size.factor, dat, mean)agg_gp$group.size.factor <-factor(agg_gp$group.size.factor)# raincloud plot:gg_all <-ggplot(dat[dat$group.size =="All sizes", ],aes(y = distance,x = type,color = type,fill = type )) +geom_signif(y_position =c(4.35, 4.75),xmin =c(1, 1),xmax =c(2, 3),annotation =c('beta:0.32~ (0.19 - 0.45)', 'beta:0.36~ (0.24 - 0.47)'), parse =TRUE,tip_length =0.02,textsize =3,size =0.5,color ="gray60",vjust =0.2 ) +# add half-violin from {ggdist} package ggdist::stat_halfeye(# fill = fill_color,alpha =0.5,# custom bandwidthadjust = .5,# adjust heightwidth = .6,.width =0,# move geom to the crightjustification =-.2,point_colour =NA ) +geom_boxplot(# fill = fill_color,width = .15,# remove outliersoutlier.shape =NA) +geom_boxplot(fill =adjustcolor("orange", alpha.f =0.1),data = agg_gp[!agg_gp$group.size %in%2:5,],aes(y = distance, x = type),color =adjustcolor("orange", alpha.f =0.6),width = .1,# remove outliersoutlier.shape =NA ) +# add justified jitter from the {gghalves} package gghalves::geom_half_point(# color = fill_color,# draw jitter on the leftside ="l",# control range of jitterrange_scale = .4,# add some transparencyalpha = .5,transformation = ggplot2::position_jitter(height =0) ) +# scale_color_viridis_d(option = "G", end = 0.8) +# scale_fill_viridis_d(option = "G",# end = 0.8,# alpha = 0.6) +scale_color_manual(values = graph_cols[2:4]) +scale_fill_manual(values = graph_cols_transp[2:4]) + gghalves::geom_half_point(data = agg_gp[agg_gp$group.size !="All sizes", ],aes(y = distance, x = type),color ="orange",# draw jitter on the leftside ="l",pch =20,size =4,# control range of jitterrange_scale = .4,# add some transparencyalpha = .5,transformation = ggplot2::position_jitter(height =0) ) +ylim(c(-0.4, 4.8)) +scale_x_discrete(labels =c("natural"="Natural","artificial"="Artificial","simulated"="Simulated" )) +theme(legend.position ="none") +labs(x ="Flight type", y ="Distance (m)") +facet_wrap(~ group.size,ncol =1,labeller =labeller(group.size =c("All sizes"="All group sizes"))) +geom_text(lineheight =0.8,data = n_grps[!n_grps$group.size.factor %in%2:6, ],aes(y = distance, x = type, label = n.labels),color ="gray60",nudge_x =0,size =3 )# add significancegg_all
Code
# raincloud plot:gg_gs <-ggplot(dat[dat$group.size %in%2:5,],aes(y = distance,x = type,color = type,fill = type )) +# add half-violin from {ggdist} package ggdist::stat_halfeye(# fill = fill_color,alpha =0.6,# custom bandwidthadjust = .5,# adjust heightwidth = .6,.width =0,# move geom to the crightjustification =-.2,point_colour =NA ) +# add justified jitter from the {gghalves} package gghalves::geom_half_point(# color = fill_color,# draw jitter on the leftside ="l",# control range of jitterrange_scale = .4,# add some transparencyalpha = .3,transformation = ggplot2::position_jitter(height =0) ) +# scale_color_viridis_d(option = "G", end = 0.8) +# scale_fill_viridis_d(option = "G",# end = 0.8,# alpha = 0.6) +scale_color_manual(values = graph_cols[2:4]) +scale_fill_manual(values = graph_cols_transp[2:4]) + gghalves::geom_half_point(data = agg_gp[agg_gp$group.size %in%2:5, ],aes(y = distance, x = type),color ="orange",# draw jitter on the leftside ="l",pch =20,size =4,# control range of jitterrange_scale = .4,# add some transparencyalpha = .5,transformation = ggplot2::position_jitter(height =0) ) +geom_boxplot(# fill = fill_color,width = .15,# remove outliersoutlier.shape =NA) +geom_boxplot(fill =adjustcolor("orange", alpha.f =0.1),data = agg_gp[agg_gp$group.size %in%2:5, ],aes(y = distance, x = type),color =adjustcolor("orange", alpha.f =0.6),width = .1,# remove outliersoutlier.shape =NA ) +ylim(c(-1.2, 4.8)) +scale_x_discrete(labels =c("natural"="Natural","artificial"="Artificial","simulated"="Simulated" )) +theme(legend.position ="none") +facet_wrap(~ group.size.factor,ncol =2,labeller =labeller(group.size.factor = custom_labels) ) +labs(x ="Flight type", y ="Pairwise flight distance (m)") +geom_text(lineheight =0.8,data = n_grps[n_grps$group.size.factor %in%2:5, ],aes(y = distance, x = type, label = n.labels),color ="gray60",nudge_x =0,size =2.3 )# gg_gspg <-plot_grid( gg_gs +theme_classic(base_size =15) +theme(legend.position ="none",axis.text.x =element_text(angle =45, hjust =1),axis.title.x =element_blank() ), gg_all +theme_classic(base_size =15) +theme(legend.position ="none", axis.title.y =element_blank()) +labs(x =""),ncol =2 )# Add a common x-axis labelpg <-ggdraw(pg) +draw_label("Group type",x =0.51,y =0.02,vjust =0,size =15 )hedFont <-"Arial"pg <- pg +theme(plot.title =element_text(size =20, family = hedFont, face ="bold"))pg
Add the group roost entry data to each simulated set and run a bayesian mixed effects model coordination on roost entry for each data set
The mean pairwise spatial distance between individuals at each frame was used as the response variable, while group type (natural, artificial or simulated) was used as predictor
Group ID and group size were included as varying effects.
The posterior distributions from the 30 models were combined into a single model fit for statistical inference
Calculate time lapse between the first individual to enter the roost for each group:
Code
group_dat <- dat[!is.infinite(dat$entry.time) &!is.na(dat$entry.time) & dat$type !="Individual", ]# get time difference between entriesgroup_dat_l <-lapply(unique(group_dat$group), function(x) {# print(x) X <- group_dat[group_dat$group == x, ] X <- X[!is.na(X$entry.time), ] X <- X[order(X$entry.time), ]# X$entry.time.diff <- X$entry.time - min(X$entry.time) X$entry.time.diff <-c(NA, X$entry.time[-1] - X$entry.time[-nrow(X)]) X <- X[!is.na(X$entry.time.diff), ]# add 1 milisecond to 0 X$entry.time.diff[X$entry.time.diff ==0] <-NA X$group.size <-if (nrow(X) >0)nrow(X) -1elsevector()return(X)})group_dat <-do.call(rbind, group_dat_l)
2.3 Simulating entry time lapse data
Simulated 30 data sets, each one with a number of simulated groups for each group size (2 to 5 individuals) similar to those in the natural group data. This approach allowed us to propagate uncertainty from random group simulations into our statistical inference.
Code
group_dat <- dat[!is.infinite(dat$entry.time) &!is.na(dat$entry.time) & dat$type !="Individual", ]# get difference to first entrygroup_dat_l <-lapply(unique(group_dat$group), function(x) {# print(x) X <- group_dat[group_dat$group == x, ] X <- X[!is.na(X$entry.time), ] X <- X[order(X$entry.time), ]# X$entry.time.diff <- X$entry.time - min(X$entry.time) X$entry.time.diff <-c(NA, X$entry.time[-1] - X$entry.time[-nrow(X)]) X <- X[!is.na(X$entry.time.diff), ] X$group.size <-if (nrow(X) >0)nrow(X) +1elsevector()return(X)})group_dat <-do.call(rbind, group_dat_l)group_dat$type <-factor(group_dat$type, levels =c("Artificial","Natural", "Simulated"))# make random groups from individual flightsindiv_dat <- dat[!is.infinite(dat$entry.time) &!is.na(dat$entry.time) & dat$type =="Individual", ]indivs <-unique(indiv_dat$Individuo)group_sizes <- group_dat$group.size[!duplicated(group_dat$Video)]# use only group sizes 2:5group_sizes <- group_sizes[group_sizes <=5]table(group_sizes)
# custom_labels <- unique(sim_groups_l[[3]]$group.size)# custom_labels <- paste("group size = ", custom_labels)# names(custom_labels) <- unique(sim_groups_l[[3]]$group.size)# table(call_rate_by_group$group.size, call_rate_by_group$type)cols <- graph_cols[2:4]cols_transp <- graph_cols_transp[2:4]names(cols) <-names(cols_transp) <-c("Natural", "Artificial", "Simulated")# raincloud plot:gg_entry <-ggplot(sim_groups_l[[1]],aes(y = entry.time.diff,x = type,color = type,fill = type )) +# add half-violin from {ggdist} package ggdist::stat_halfeye(# fill = fill_color,alpha =0.5,# custom bandwidthadjust = .5,# adjust heightwidth = .6,.width =0,# move geom to the crightjustification =-.2,point_colour =NA ) +geom_boxplot(# fill = fill_color,width = .15,# remove outliersoutlier.shape =NA) +# add justified jitter from the {gghalves} package gghalves::geom_half_point(# color = fill_color,# draw jitter on the leftside ="l",# control range of jitterrange_scale = .4,# add some transparencyalpha = .5,transformation = ggplot2::position_jitter(height =0) ) +scale_color_manual(values = cols) +scale_fill_manual(values = cols_transp) +# scale_color_viridis_d(option = "G", end = 0.8) +# scale_fill_viridis_d(option = "G",# end = 0.8,# alpha = 0.6) +# ylim(c(-0.1, 4.8)) +scale_x_discrete(labels =c("natural"="Natural","artificial"="Artificial","simulated"="Simulated" )) +theme(legend.position ="none") +labs(x ="Group type", y ="Time difference (s)")agg_grp_time <-aggregate(entry.time.diff ~ type + group, sim_groups_l[[1]], mean)n_grps <-aggregate(group ~ type + group.size, sim_groups_l[[1]], function(x)length(unique(x)))n_grps$indivs <- n_grps$group * n_grps$group.sizeagg_grps <-aggregate(group ~ type, n_grps, sum)agg_grps$indivs <-aggregate(indivs ~ type, n_grps, sum)$indivsagg_grps$n.labels <-paste0(agg_grps$group, " / ", agg_grps$indivs)agg_grps$n.labels[1] <-" n = 13 groups /\n 50 individuals"agg_grps$distance <--15gg_entry <- gg_entry +ylim(c(-20, 210)) +geom_signif(y_position =c(170, 190, 210),xmin =c(1, 1, 2),xmax =c(3, 2, 3),annotation =c('beta:-1.71~ (-2.31 - -1.103)', 'beta:-0.99~ (-1.59 - -0.40)', 'beta:-0.72 ~ (-1.30 - -0.12)'), parse =TRUE,tip_length =0.02,textsize =3,size =0.5,color ="gray60",vjust =0.2 ) + gghalves::geom_half_point(data = agg_grp_time,aes(y = entry.time.diff, x = type),color ="orange",# draw jitter on the leftside ="l",pch =20,size =4,# control range of jitterrange_scale = .4,# add some transparencyalpha = .5,transformation = ggplot2::position_jitter(height =0) ) +geom_boxplot(fill =adjustcolor("orange", alpha.f =0.1),data = agg_grp_time,aes(y = entry.time.diff, x = type),color =adjustcolor("orange", alpha.f =0.6),width = .1,# remove outliersoutlier.shape =NA ) +geom_text(lineheight =0.8,data = agg_grps,aes(y = distance, x = type, label = n.labels),color ="gray60",nudge_x =0,size =4 ) +theme_classic(base_size =12) +theme(legend.position ="none")ggsave("./output/roost_entry_coordination_1_panel.png", gg_entry, grDevices::png,width =4.5,height =3.5,dpi =300)
By group size:
Code
# raincloud plot:ggplot(sim_groups_l[[1]],aes(y = entry.time.diff,x = type,color = type,fill = type )) +# add half-violin from {ggdist} package ggdist::stat_halfeye(# fill = fill_color,alpha =0.5,# custom bandwidthadjust = .5,# adjust heightwidth = .6,.width =0,# move geom to the crightjustification =-.2,point_colour =NA ) +geom_boxplot(# fill = fill_color,width = .15,# remove outliersoutlier.shape =NA) +# add justified jitter from the {gghalves} package gghalves::geom_half_point(# color = fill_color,# draw jitter on the leftside ="l",# control range of jitterrange_scale = .4,# add some transparencyalpha = .5,transformation = ggplot2::position_jitter(height =0) ) +scale_color_viridis_d(option ="G", end =0.8) +scale_fill_viridis_d(option ="G",end =0.8,alpha =0.6) +# ylim(c(-0.1, 4.8)) +scale_x_discrete(labels =c("natural"="Natural","artificial"="Artificial","simulated"="Simulated" )) +theme(legend.position ="none") +facet_wrap(~ group.size.factor,ncol =2,labeller =labeller(group.size.factor = custom_labels) ) +labs(x ="Flight type", y ="Time difference (s)")
2.5 Regression model
Add the group roost entry data to each simulated set and run a bayesian mixed effects model coordination on roost entry for each data set
The difference in time between when entering the roost between consecutive individuals was used as the response variable (modeled with a normal distribution), while group type (natural, artificial or simulated) was used as predictor
Group ID was included as varying effect
The posterior distributions from the 30 models were combined into a single model fit for statistical inference
mod <-brm(formula = calls |resp_rate(flight.time) ~ type +mo(group.size.f) + (1| indiv), iter = iter, thin =1, data = call_rate_indiv, family =negbinomial(),silent =2, chains = chains, cores = chains, control =list(adapt_delta =0.99,max_treedepth =15), file ="./data/processed/individual_call_rate_solo_vs_group_group_size_monot_negbinomial",file_refit ="always")mod <-add_criterion(mod, criterion ="loo")null_mod <-brm(formula = calls |resp_rate(flight.time) ~1+ (1| indiv), iter = iter, thin =1, data = call_rate_indiv, family =negbinomial(),silent =2, chains = chains, cores = chains, control =list(adapt_delta =0.99,max_treedepth =15), file ="./data/processed/null_model_individual_call_rate_solo_vs_group_group_size_monot_negbinomial",file_refit ="always")null_mod <-add_criterion(null_mod, criterion ="loo")
3.2.2.2.1 Results
3.2.2.2.1.1 Model performance vs null model
Code
mod <-readRDS("./data/processed/individual_call_rate_solo_vs_group_group_size_monot.rds")null_mod <-readRDS("./data/processed/null_model_individual_call_rate_solo_vs_group_group_size_monot.rds")loo_compare(mod, null_mod)
mod <-brm(formula =log(cv.gaps +1) ~ type +mo(group.size.f) + (1| group), iter = iter, thin =1, data = cv_gaps, family =gaussian(),silent =2, chains = chains, cores = chains, file ="./data/processed/cv_gaps_solo_vs_group",file_refit ="always")mod <-add_criterion(mod, criterion ="loo")null_mod <-brm(formula =log(cv.gaps +1) ~1+ (1| group), iter = iter,thin =1, data = cv_gaps, family =gaussian(), silent =2, chains = chains,cores = chains, file ="./data/processed/null_cv_gaps_solo_vs_group",file_refit ="always")null_mod <-add_criterion(null_mod, criterion ="loo")beepr::beep(2)
3.4.1.2.1 Results
3.4.1.2.1.1 Model performance vs null model
Code
mod <-readRDS("./data/processed/cv_gaps_solo_vs_group.rds")null_mod <-readRDS("./data/processed/null_cv_gaps_solo_vs_group.rds")loo_compare(mod, null_mod)
Response modeled with an exponential distribution.
Code
mod <-brm(formula =log(gaps +1) ~ type +mo(group.size.f) + (1| indiv) + (1| group), iter = iter, thin =1, data = sub_gaps_indiv,family =gaussian(), silent =2, chains = chains, cores = chains,control =list(adapt_delta =0.99, max_treedepth =15), file_refit ="always",file ="./data/processed/group_gaps_solo_vs_group_group_size_monot")custom_ppc(mod, group ="type")mod <-add_criterion(mod, criterion ="loo")null_mod <-brm(formula =log(gaps +1) ~1+ (1| indiv) + (1| group),iter = iter, thin =1, data = sub_gaps_indiv, family =gaussian(),silent =2, chains = chains, cores = chains, control =list(adapt_delta =0.99,max_treedepth =15), file_refit ="always", file ="./data/processed/null_group_gaps_solo_vs_group_group_size_monot")null_mod <-add_criterion(null_mod, criterion ="loo")beepr::beep(2)
3.4.2.2.1 Results
3.4.2.2.2 Model performance vs null model
Code
mod <-readRDS("./data/processed/group_gaps_solo_vs_group_group_size_monot.rds")null_mod <-readRDS("./data/processed/null_group_gaps_solo_vs_group_group_size_monot.rds")comp <-loo_compare(mod, null_mod)print(comp)
elpd_diff se_diff
mod 0.0 0.0
null_mod -165.3 24.3
mod <-brm(formula = norm.markov.entropy ~ type, iter = iter, thin =1,data = seq_entropy_data, family =zero_one_inflated_beta(), silent =2,chains = chains, cores = chains, control =list(adapt_delta =0.99,max_treedepth =15), file_refit ="always", file ="./data/processed/entropy_by_experiment_type_group_size_monot")custom_ppc(mod, group ="type")mod <-add_criterion(mod, criterion ="loo")null_mod <-brm(formula = norm.markov.entropy ~1, iter = iter, thin =1,data = seq_entropy_data, family =zero_one_inflated_beta(), silent =2,chains = chains, cores = chains, control =list(adapt_delta =0.99,max_treedepth =15), file_refit ="always", file ="./data/processed/null_model_entropy_by_experiment_type_group_size_monot")null_mod <-add_criterion(null_mod, criterion ="loo")
3.5.2.1 Results
3.5.2.1.1 Model performance vs null model
Code
mod <-readRDS("./data/processed/entropy_by_experiment_type_group_size_monot.rds")null_mod <-readRDS("./data/processed/null_model_entropy_by_experiment_type_group_size_monot.rds")comp <-loo_compare(mod, null_mod)print(comp)
No difference in predictability between natural and artificial groups
3.6 Number of calling individuals per experiment
3.6.1 Plot raw data
All groups combined:
Code
# raincloud plot:ggplot(seq_entropy_data,aes(y = prop.callers,x = type,color = type,fill = type )) +# add half-violin from {ggdist} package ggdist::stat_halfeye(# fill = fill_color,alpha =0.5,# custom bandwidthadjust = .5,# adjust heightwidth = .6,.width =0,# move geom to the crightjustification =-.2,point_colour =NA ) +geom_boxplot(# fill = fill_color,width = .15,# remove outliersoutlier.shape =NA) +# add justified jitter from the {gghalves} package gghalves::geom_half_point(# color = fill_color,# draw jitter on the leftside ="l",# control range of jitterrange_scale = .4,# add some transparencyalpha = .5,transformation = ggplot2::position_jitter(height =0) ) +scale_color_viridis_d(option ="G", end =0.8) +scale_fill_viridis_d(option ="G",end =0.8,alpha =0.6) +# ylim(c(-0.1, 4.8)) +# scale_x_discrete(labels = c(# "natural" = "Natural",# "artificial" = "Artificial",# "simulated" = "Simulated"# )) +theme(legend.position ="none") +labs(x ="Flight type", y ="Proportion of unique callers")
Warning in bandwidth_dpi(): Bandwidth calculation failed.
→ Falling back to `bandwidth_nrd0()`.
ℹ This often occurs when a sample contains many duplicates, which suggests that
a dotplot (e.g., `geom_dots()`) or histogram (e.g., `density_histogram()`,
`stat_slab(density = 'histogram')`, or `stat_histinterval()`) may better
represent the data.
Caused by error in `bw.SJ()`:
! sample is too sparse to find TD
By group size:
Code
# raincloud plot:ggplot(seq_entropy_data[seq_entropy_data$group.size <6, ],aes(y = prop.callers,x = type,color = type,fill = type )) +# add half-violin from {ggdist} package ggdist::stat_halfeye(# fill = fill_color,alpha =0.5,# custom bandwidthadjust = .5,# adjust heightwidth = .6,.width =0,# move geom to the crightjustification =-.2,point_colour =NA ) +geom_boxplot(# fill = fill_color,width = .15,# remove outliersoutlier.shape =NA) +# add justified jitter from the {gghalves} package gghalves::geom_half_point(# color = fill_color,# draw jitter on the leftside ="l",# control range of jitterrange_scale = .4,# add some transparencyalpha = .5,transformation = ggplot2::position_jitter(height =0) ) +scale_color_viridis_d(option ="G", end =0.8) +scale_fill_viridis_d(option ="G",end =0.8,alpha =0.6) +# ylim(c(-0.1, 4.8)) +# scale_x_discrete(labels = c(# "natural" = "Natural",# "artificial" = "Artificial",# "simulated" = "Simulated"# )) +theme(legend.position ="none") +facet_wrap(~ group.size,ncol =2 ) +labs(x ="Flight type", y ="Proporiton of unique callers")
3.6.2 Regression model
Code
mod <-brm(formula = prop.callers ~ type + (1| group.size.f), iter = iter,thin =1, data = seq_entropy_data, family =zero_one_inflated_beta(),silent =2, chains = chains, cores = chains, control =list(adapt_delta =0.99,max_treedepth =15), file_refit ="always", file ="./data/processed/unique.callers_by_experiment_type_group_size_monot")custom_ppc(mod, group ="type")mod <-add_criterion(mod, criterion ="loo")null_mod <-brm(formula = prop.callers ~1+ (1| group.size.f), iter = iter,thin =1, data = seq_entropy_data, family =zero_one_inflated_beta(),silent =2, chains = chains, cores = chains, control =list(adapt_delta =0.99,max_treedepth =15), file_refit ="always", file ="./data/processed/null_model_unique.callers_by_experiment_type_group_size_monot")null_mod <-add_criterion(null_mod, criterion ="loo")
3.6.2.1 Results
3.6.2.1.1 Model performance vs null model
Code
mod <-readRDS("./data/processed/unique.callers_by_experiment_type_group_size_monot.rds")null_mod <-readRDS("./data/processed/null_model_unique.callers_by_experiment_type_group_size_monot.rds")comp <-loo_compare(mod, null_mod)print(comp)
