opts_knit$set(root.dir = "..")

opts_chunk$set(fig.width = 12, fig.height = 8)

theme_set(theme_classic(base_size = 34))

# print results brms models
summary_brm_model <- function(x, gsub.pattern = NULL, gsub.replacement = NULL, cols = viridis(10, alpha = 0.7)){
  
  summ <- summary(x)$fixed
  fit <- x$fit  
  betas <- grep("^b_", names(fit@sim$samples[[1]]), value = TRUE)  
  hdis <- t(sapply(betas, function(y)   hdi(fit@sim$samples[[1]][[y]]))
)
  summ_table <- data.frame(summ, hdis, iterations = attr(fit, "stan_args")[[1]]$iter, chains = length(attr(fit, "stan_args")))
  summ_table <- summ_table[rownames(summ_table) != "Intercept", c("Estimate", "Rhat", "Bulk_ESS", "l.95..CI", "u.95..CI", "iterations", "chains")]
  
  summ_table <- as.data.frame(summ_table)  
  summ_table$Rhat <- round(summ_table$Rhat, digits = 3)  
  summ_table$CI_low <- round(unlist(summ_table$l.95..CI), digits = 3)  
  summ_table$CI_high <- round(unlist(summ_table$u.95..CI), digits = 3)  
  summ_table$l.95..CI <- summ_table$u.95..CI <- NULL
  
    out <- lapply(betas, function(y)  data.frame(variable = y, value = sort(fit@sim$samples[[1]][[y]], decreasing = FALSE)))
  
  posteriors <- do.call(rbind, out)
  posteriors <- posteriors[posteriors$variable != "b_Intercept", ]
  
  if (!is.null(gsub.pattern) & !is.null(gsub.replacement))
    posteriors$variable <- gsub(pattern = gsub.pattern, replacement = gsub.replacement, posteriors$variable)
  
   gg <- ggplot(data = posteriors, aes(y = variable, x = value, fill = stat(quantile))) + 
  geom_density_ridges_gradient(quantile_lines = TRUE, quantile_fun = HDInterval::hdi,  vline_linetype = 2) + 
  theme_classic() + 
     xlim(range(c(min(summ_table[ , "CI_low"]), max(summ_table[ , "CI_high"])), 0)) +
  geom_vline(xintercept = 0, col = "red") + 
  scale_fill_manual(values = c("transparent", "lightblue", "transparent"), guide = "none") + ggtitle(x$formula)

  if (!is.null(gsub.pattern) & !is.null(gsub.replacement))
    rownames(summ_table) <- gsub(pattern = gsub.pattern, replacement = gsub.replacement, rownames(summ_table))
   
  summ_table$Rhat <- ifelse(summ_table$Rhat > 1.05, cell_spec(summ_table$Rhat, "html", color ="white", background = "red", bold = TRUE,  font_size = 12),  cell_spec(summ_table$Rhat, "html"))
  
  signif <- summ_table[,"CI_low"] * summ_table[,"CI_high"] > 0
  
  df1 <- kbl(summ_table, row.names = TRUE, escape = FALSE, format = "html", digits = 3)
    
  df1 <- row_spec(kable_input = df1,row =  which(summ_table$CI_low * summ_table$CI_high > 0), background = adjustcolor(cols[9], alpha.f = 0.3))

  df1 <- kable_styling(df1, bootstrap_options = c("striped", "hover", "condensed", "responsive"), full_width = FALSE, font_size = 12)
  
  print(df1)
  
  print(gg)

  }

Data description

dat <- as.data.frame(read_excel(path = "./data/raw/Anexo 1_Entrada a refugios-2022.xlsx"))

dat$entry.time <- as.numeric(dat$`Tiempo real`)

## Warning: NAs introduzidos por coerção

# dat$entry.time
# 
# str(dat)

cat("Number of individuals per type:")

## Number of individuals per type:

(table(dat$`Tipo grupo`))

## 
## Individual      Mixto       Real 
##         35        118        159

cat("Number of individuals per type removing missing data (NAs and Inf:")

## Number of individuals per type removing missing data (NAs and Inf:

(table(dat$`Tipo grupo`[!is.infinite(dat$entry.time) & !is.na(dat$entry.time)]))

## 
## Individual      Mixto       Real 
##         17         65         95

Stats

Entry time difference

group_dat <- dat[!is.infinite(dat$entry.time) & !is.na(dat$entry.time) & dat$`Tipo grupo` != "Individual", ]

# get difference to first entry
group_dat_l <- lapply(unique(group_dat$Grupo), function(x){
  # print(x)
  X <- group_dat[group_dat$Grupo == x, ]
  X$entry.time.diff <- X$entry.time - min(X$entry.time)
  X <- X[-which.min(X$entry.time.diff), ]
  
  X$group.size <- if (nrow(X) > 0) nrow(X) - 1 else vector()
  return(X)
  })

group_dat <- do.call(rbind, group_dat_l)

group_dat$Group <- group_dat$Grupo
group_dat$entry.time.diff[group_dat$entry.time.diff == 0] <- 0.001
group_dat$type <- factor(ifelse(group_dat$`Tipo grupo` == "Mixto", "Aritficial", "Real"), levels = c("Aritficial", "Real"))

# individual plot solo vs individual in group
ggplot(group_dat, aes(x = type, y = log(entry.time.diff))) + 
  geom_boxplot() +
  # geom_point(size = 2) +
  # scale_color_viridis_c(alpha = 0.8) +
  labs(x = "Group type", y = "log(time after first entry (s))") + 
  # facet_grid( ~ experiment_f) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

chains <- 1
iter <- 5000

model_formulas <- c("entry.time.diff ~ type + (1 | Group) + (1 | group.size)", "entry.time.diff ~ 1 + (1 | Group) + (1 | group.size)")

brms_models <- lapply(model_formulas, function(x){

  mod <- brm(
          formula = x,
          iter = iter,
          thin = 1,
          data = group_dat,
          family = lognormal(),
          silent = 2,
          chains = chains,
          cores = chains
          )
  
  mod <- add_criterion(mod, c("loo"))

  return(mod)
  })

names(brms_models) <- model_formulas

saveRDS(brms_models, "./data/processed/roost_entry_coordination_models.RDS")

brms_models <- readRDS("./data/processed/roost_entry_coordination_models.RDS")

cat(paste(length(brms_models),"models evaluated:\n"))

2 models evaluated:

for(i in names(brms_models))
  cat(paste("-  ", i, "\n"))

entry.time.diff ~ type + (1 | Group) + (1 | group.size)
entry.time.diff ~ 1 + (1 | Group) + (1 | group.size)

comp_mods <- loo_compare(brms_models[[1]], brms_models[[2]], model_names = names(brms_models))


df1 <- kbl(comp_mods, row.names = TRUE, escape = FALSE, format = "html", digits = 3)

cat("Compare models:")

Compare models:

 kable_styling(df1, bootstrap_options = c("striped", "hover", "condensed", "responsive"), full_width = FALSE, font_size = 12)

	elpd_diff	se_diff	elpd_loo	se_elpd_loo	p_loo	se_p_loo	looic	se_looic
entry.time.diff ~ 1 + (1 \| Group) + (1 \| group.size)	0.000	0.000	-453.128	27.109	20.834	3.708	906.256	54.217
entry.time.diff ~ type + (1 \| Group) + (1 \| group.size)	-0.781	2.253	-453.909	26.927	18.435	3.646	907.818	53.853

cat("Best model:\n")

Best model:

cat(paste("-  ", rownames(comp_mods)[1], "\n"))

entry.time.diff ~ 1 + (1 | Group) + (1 | group.size)

# best model
if (!grepl("1 +", rownames(comp_mods)[1], fixed = TRUE))
summary_brm_model(brms_models[[rownames(comp_mods)[1]]], gsub.pattern = "experiment.type", gsub.replacement = "solo_vs_")

Takeaways

No differences in roost entry coordination between real and artificial groups

Simulation entry time difference

group_dat <- dat[!is.infinite(dat$entry.time) & !is.na(dat$entry.time) & dat$`Tipo grupo` != "Individual", ]

# get difference to first entry
group_dat_l <- lapply(unique(group_dat$Grupo), function(x){
  # print(x)
  X <- group_dat[group_dat$Grupo == x, ]
  X$entry.time.diff <- X$entry.time - min(X$entry.time)
  X <- X[-which.min(X$entry.time.diff), ]
  
  X$group.size <- if (nrow(X) > 0) nrow(X) + 1 else vector()
  return(X)
  })

group_dat <- do.call(rbind, group_dat_l)



group_dat$Group <- group_dat$Grupo
group_dat$entry.time.diff[group_dat$entry.time.diff == 0] <- 0.001
group_dat$type <- factor(ifelse(group_dat$`Tipo grupo` == "Mixto", "Artificial", "Real"), levels = c("Simulated", "Artificial", "Real"))

# make random groups from individual flights
indiv_dat <- dat[!is.infinite(dat$entry.time) & !is.na(dat$entry.time) & dat$`Tipo grupo` == "Individual", ]

indivs <- unique(indiv_dat$Individuo)


sim_groups_l <- pblapply(1:1000, function(x){
  
 g_size <- sample(group_dat$group.size)
  
 n_group <- sum(cumsum(g_size) <= length(indivs))  
 
 g_size <- g_size[1:n_group]
 
 sampled_indivs <- sample(indivs, sum(g_size))
 
 indivs_split <- split(sampled_indivs, f = unlist(sapply(1:n_group, function(x) rep(x, g_size[x]))))
 
 sub_sim_groups_l <- lapply(1:length(indivs_split), function(y) {
   
   sim_group <- indiv_dat[indiv_dat$Individuo %in% indivs_split[[y]], ]
   sim_group$Group <- paste("sim",x, y, sep = "-")
 
   sim_group$entry.time.diff <- sim_group$entry.time - min(sim_group$entry.time)
   sim_group <- sim_group[-which.min(sim_group$entry.time.diff), ]
   sim_group$group.size <- nrow(sim_group) + 1
   sim_group$type <- factor("Simulated")
   
   out <- rbind(sim_group, group_dat)
   
   return(out)
   }
  )
 
 sub_sim_group <- do.call(rbind, sub_sim_groups_l)
 
 return(sub_sim_group)
})


sim_group_n <- sapply(sim_groups_l, function(x) length(unique(x$Group[x$type == "Simulated"])))

table(sim_group_n)

# keep those with same number of simulated groups
sim_groups_l <- sim_groups_l[which(sim_group_n == 5)]

sim_groups_l <- lapply(sim_groups_l, function(x){
  
  Y <- x[x$type == "Simulated", ]
  Y$Group <- as.numeric(as.factor(Y$Group)) 
  Y <- rbind(Y, x[x$type != "Simulated", ])
  
  return(Y)
  })

# keep only 30
sim_groups_l <- sim_groups_l[1:30]

sim_group_dat <- do.call(rbind, lapply(sim_groups_l, function(x) x[x$type == "Simulated", ]))

sim_group_dat <- rbind(sim_group_dat, group_dat)


# individual plot solo vs individual in group
ggplot(sim_group_dat, aes(x = type, y = log(entry.time.diff))) + 
  geom_boxplot() +
  # geom_point(size = 2) +
  # scale_color_viridis_c(alpha = 0.8) +
  labs(x = "Flight", y = "log(time after first entry (s))") + 
  # facet_grid( ~ experiment_f) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

saveRDS(sim_groups_l, "./data/processed/simulated_data_entry_time.RDS")

sim_groups_l <- readRDS("./data/processed/simulated_data_entry_time.RDS")

chains <- 1
iter <- 10000

mods <- brm_multiple(
          formula = entry.time.diff ~ type + (1 | Group) + (1 | group.size),
          iter = iter,
          thin = 1,
          data = sim_groups_l,
          family = lognormal(),
          silent = 2,
          chains = chains,
          cores = chains,
          combine = FALSE
          )

combined_mod <- combine_models(mlist = mods, check_data = FALSE)


saveRDS(combined_mod, "./data/processed/simulated_data_entry_time_models.RDS")

combined_mod <- readRDS("./data/processed/simulated_data_entry_time_models.RDS")

summary_brm_model(combined_mod)

## Warning: There were 816 divergent transitions after warmup. Increasing
## adapt_delta above 0.8 may help. See http://mc-stan.org/misc/
## warnings.html#divergent-transitions-after-warmup

	Estimate	Rhat	Bulk_ESS	iterations	chains	CI_low	CI_high
typeArtificial	-3.221	1.018	28996.700	10000	30	-5.884	-0.609
typeReal	-4.524	1.027	1265.216	10000	30	-6.752	-2.251

## Picking joint bandwidth of 0.166

Real vs artificial:

contrasts <- c(real_vs_artificial = "typeReal - typeArtificial = 0")
hypothesis(combined_mod, contrasts)

## Hypothesis Tests for class b:
##           Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio Post.Prob
## 1 real_vs_artificial     -1.3      0.92    -3.09     0.49         NA        NA
##   Star
## 1     
## ---
## 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 95%;
## for two-sided hypotheses, the value tested against lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.

Takeaways

Both real and artificial groups enter more sync’ed than simulated groups
No differences in roost entry coordination between real and artificial groups

Propensity to enter the roost

group_dat <- dat[dat$`Tipo grupo` != "Individual", ]
group_dat$indiv.entry <- ifelse(is.na(group_dat$entry.time) | is.infinite(group_dat$entry.time), 0, 1)


group_dat$Group <- group_dat$Grupo
group_dat$type <- factor(ifelse(group_dat$`Tipo grupo` == "Mixto", "Aritficial", "Real"), levels = c("Aritficial", "Real"))  

agg_dat <- aggregate(indiv.entry ~ Group + type, group_dat, FUN = sum)
agg_dat$group.size <- aggregate(indiv.entry ~ Group + type, group_dat, FUN = length)$indiv.entry


agg_dat$entry <- "Partial"
agg_dat$entry[agg_dat$group.size == agg_dat$indiv.entry] <- "Full"
agg_dat$entry[agg_dat$indiv.entry == 0] <- "None"

count_dat <- aggregate(Group ~ entry + type, agg_dat, length)
count_dat$entry <- factor(count_dat$entry, levels = c("None", "Partial", "Full"))

# Change barplot fill colors by groups
ggplot(count_dat, aes(x=entry, y=Group, fill=type)) +
  geom_bar(stat="identity", position=position_dodge()) + 
  scale_fill_viridis_d(begin = 0.3, end = 0.8) +
  theme_classic()

# mean centering
agg_dat$group.size <- agg_dat$group.size - mean(agg_dat$group.size)

chains <- 1
iter <- 10000

agg_dat$entry[agg_dat$entry == "None"] <-  "1_None"

model_formulas <- c("entry ~ type + (1 | Group) + ( 1 | group.size)")

brms_models <- lapply(model_formulas, function(x){
print(x)
  mod <- brm(data = agg_dat, 
            family = categorical, 
            formula = x,
            iter = iter,
            chains = chains, 
            cores = chains
              )
  
  try(mod <- add_criterion(mod, c("loo")))

  return(mod)
  })

names(brms_models) <- model_formulas

saveRDS(brms_models, "./data/processed/roost_entry_models.RDS")

brms_models <- readRDS("./data/processed/roost_entry_models.RDS")

summary_brm_model(brms_models$`entry ~ type + (1 | Group) + ( 1 | group.size)`)

## Warning: There were 2 divergent transitions after warmup. Increasing adapt_delta
## above 0.8 may help. See http://mc-stan.org/misc/warnings.html#divergent-
## transitions-after-warmup

	Estimate	Rhat	Bulk_ESS	iterations	chains	CI_low	CI_high
muFull_Intercept	-0.693	1	560.074	10000	1	-16.059	5.511
muPartial_Intercept	6.370	1	1342.146	10000	1	1.466	17.166
muFull_typeReal	2.479	1.001	577.631	10000	1	-5.492	27.256
muPartial_typeReal	-12.342	1	1103.571	10000	1	-32.277	-3.617

## Picking joint bandwidth of 0.525

contrasts <- c(one = "muFull_Intercept - muFull_typeReal = 0")
hypothesis(brms_models$`entry ~ type + (1 | Group) + ( 1 | group.size)`, contrasts)

## Hypothesis Tests for class b:
##   Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio Post.Prob Star
## 1        one    -3.17     12.28   -42.72    10.63         NA        NA     
## ---
## 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 95%;
## for two-sided hypotheses, the value tested against lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.

brms_models <- readRDS("./data/processed/roost_entry_models.RDS")

cat(paste(length(brms_models),"models evaluated:\n"))

for(i in names(brms_models))
  cat(paste("-  ", i, "\n"))

comp_mods <- loo_compare(brms_models[[1]], brms_models[[2]],  brms_models[[3]], model_names = names(brms_models))


df1 <- kbl(comp_mods, row.names = TRUE, escape = FALSE, format = "html", digits = 3)

cat("Compare models:")
  
 kable_styling(df1, bootstrap_options = c("striped", "hover", "condensed", "responsive"), full_width = FALSE, font_size = 12)

cat("Best model:\n")
cat(paste("-  ", rownames(comp_mods)[1], "\n"))

# best model
if (!grepl("1 +", rownames(comp_mods)[1], fixed = TRUE))
  summary_brm_model(brms_models[[rownames(comp_mods)[1]]], gsub.pattern = "experiment.type", gsub.replacement = "solo_vs_")

Full entry of real vs artificial groups:

real_distances_l <- readRDS("./data/processed/real_trajectory_mean_distances.RDS")

Takeaways

No clear differences between real and artificial groups in propensity to enter roost (full group entry)

Session information

## R version 4.1.0 (2021-05-18)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.2 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/atlas/libblas.so.3.10.3
## LAPACK: /usr/lib/x86_64-linux-gnu/atlas/liblapack.so.3.10.3
## 
## locale:
##  [1] LC_CTYPE=pt_BR.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=es_CR.UTF-8        LC_COLLATE=pt_BR.UTF-8    
##  [5] LC_MONETARY=es_CR.UTF-8    LC_MESSAGES=pt_BR.UTF-8   
##  [7] LC_PAPER=es_CR.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=es_CR.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] tidybayes_3.0.1   ggridges_0.5.3    kableExtra_1.3.1  brms_2.16.3      
##  [5] Rcpp_1.0.8.3      pbapply_1.5-0     readxl_1.3.1      knitr_1.37       
##  [9] ggplot2_3.3.5     viridis_0.6.2     viridisLite_0.4.0
## 
## loaded via a namespace (and not attached):
##   [1] minqa_1.2.4          colorspace_2.0-2     ellipsis_0.3.2      
##   [4] rsconnect_0.8.18     markdown_1.1         base64enc_0.1-3     
##   [7] rstudioapi_0.13      farver_2.1.0         rstan_2.21.2        
##  [10] svUnit_1.0.6         DT_0.18              fansi_1.0.0         
##  [13] mvtnorm_1.1-2        diffobj_0.3.4        xml2_1.3.2          
##  [16] bridgesampling_1.1-2 codetools_0.2-18     splines_4.1.0       
##  [19] shinythemes_1.2.0    bayesplot_1.8.1      projpred_2.0.2      
##  [22] jsonlite_1.7.2       nloptr_1.2.2.2       ggdist_3.0.0        
##  [25] shiny_1.6.0          httr_1.4.2           compiler_4.1.0      
##  [28] backports_1.3.0      assertthat_0.2.1     Matrix_1.3-4        
##  [31] fastmap_1.1.0        cli_3.1.0            later_1.2.0         
##  [34] htmltools_0.5.2      prettyunits_1.1.1    tools_4.1.0         
##  [37] igraph_1.2.6         coda_0.19-4          gtable_0.3.0        
##  [40] glue_1.6.2           reshape2_1.4.4       dplyr_1.0.7         
##  [43] posterior_1.1.0      V8_3.4.2             cellranger_1.1.0    
##  [46] jquerylib_0.1.4      vctrs_0.3.8          nlme_3.1-152        
##  [49] crosstalk_1.1.1      tensorA_0.36.2       xfun_0.30           
##  [52] stringr_1.4.0        ps_1.6.0             rvest_1.0.0         
##  [55] lme4_1.1-27.1        mime_0.11            miniUI_0.1.1.1      
##  [58] lifecycle_1.0.1      gtools_3.9.2         MASS_7.3-54         
##  [61] zoo_1.8-9            scales_1.1.1         colourpicker_1.1.0  
##  [64] promises_1.2.0.1     Brobdingnag_1.2-6    parallel_4.1.0      
##  [67] inline_0.3.19        shinystan_2.5.0      gamm4_0.2-6         
##  [70] yaml_2.3.5           curl_4.3.2           gridExtra_2.3       
##  [73] loo_2.4.1.9000       StanHeaders_2.21.0-7 sass_0.4.0          
##  [76] stringi_1.7.6        highr_0.9            dygraphs_1.1.1.6    
##  [79] checkmate_2.0.0      boot_1.3-28          pkgbuild_1.2.0      
##  [82] rlang_1.0.1          pkgconfig_2.0.3      matrixStats_0.61.0  
##  [85] HDInterval_0.2.2     distributional_0.2.2 evaluate_0.15       
##  [88] lattice_0.20-44      purrr_0.3.4          labeling_0.4.2      
##  [91] rstantools_2.1.1     htmlwidgets_1.5.3    tidyselect_1.1.1    
##  [94] processx_3.5.2       plyr_1.8.6           magrittr_2.0.2      
##  [97] R6_2.5.1             generics_0.1.0       DBI_1.1.1           
## [100] mgcv_1.8-36          pillar_1.6.4         withr_2.4.3         
## [103] xts_0.12.1           abind_1.4-5          tibble_3.1.6        
## [106] crayon_1.5.0         arrayhelpers_1.1-0   utf8_1.2.2          
## [109] rmarkdown_2.9        grid_4.1.0           callr_3.7.0         
## [112] threejs_0.3.3        webshot_0.5.2        digest_0.6.29       
## [115] xtable_1.8-4         tidyr_1.1.3          httpuv_1.6.1        
## [118] RcppParallel_5.1.4   stats4_4.1.0         munsell_0.5.0       
## [121] bslib_0.2.5.1        shinyjs_2.0.0

Roost entry coordination

Group flight coordination in Thyroptera

Marcelo Araya-Salas, PhD & Gloriana Chaverri, PhD

06-04-2022

Data description

Stats

Entry time difference

Takeaways

Simulation entry time difference

Takeaways

Propensity to enter the roost

Takeaways