Analysis of Experiment 1

Table format

outcome_tally %>% 
  kable(digits = 3) %>% kable_styling() %>%
  scroll_box(height = "380px")

Plotted as a bar chart

all_hives <- outcome_tally %>%
  group_by(treatment, outcome) %>%
  summarise(n = sum(n)) %>%
  ungroup() %>% mutate(hive = "All hives")

outcome_tally %>%
  group_by(treatment, hive, outcome) %>%
  summarise(n = sum(n)) %>%
  ungroup() %>%
  mutate(hive = paste(hive, "hive")) %>%
  full_join(all_hives, by = c("treatment", "hive", "outcome", "n")) %>%
  mutate(prop = n / sum(n)) %>%
  ggplot(aes(outcome, n , fill = treatment)) + 
  geom_bar(stat = "identity", position = "dodge") + 
  facet_wrap( ~ hive, nrow = 1, scales = "free_y") + 
  scale_fill_brewer(palette = "Set2", name = "Treatment") + 
  xlab("Outcome") + ylab("Number of bees") + 
  theme(axis.text.x = element_text(angle = 45, hjust=1))

Run the models

Fit a multinomial logisitic model, with 3 possible outcomes describing what happened to each bee introduced to the hive: stayed inside, left of its own volition, or forced out by the other workers. To assess the effects of our predictor variables, we compare 5 models with different fixed factors, ranking them by posterior model probability.

if(!file.exists("output/model_probabilities_expt2.rds")){
  my_prior <- c(set_prior("normal(0, 2)", class = "b"))
  
  run_brms_categorical <- function(formula){
    formula <- gsub(" ", "", paste("outcome_numeric ~", formula))
    brm(formula, 
        data = data_for_categorical_model, 
        prior = my_prior,
        family = "categorical", 
        save_all_pars = TRUE,
        control = list(adapt_delta = 0.999, max_treedepth = 15),
        chains = 4, cores = 4, iter = 40000, seed = 1)
  }
  
  formula_list <- c(
    "treatment:hive + treatment + hive + observation_time_minutes",
    "treatment + hive + observation_time_minutes",
    "treatment + observation_time_minutes",
    "hive + observation_time_minutes",
    "observation_time_minutes")
  
  # run all 5 models, sequentially
  model_list_expt2 <- lapply(formula_list, run_brms_categorical)
  
  # Estimate the log marginal likelihood by bridge sampling, and use to determine the posterior model probabilities
  # See: vignette("bridgesampling_example_stan")
  model_probabilities <- post_prob(model_list_expt2[[1]], 
                                   model_list_expt2[[2]], 
                                   model_list_expt2[[3]], 
                                   model_list_expt2[[4]], 
                                   model_list_expt2[[5]])
  top_index <- which.max(unname(model_probabilities))
  
  model_probabilities <- tibble(
    `Fixed effects` = formula_list, 
    `Posterior model probability` = unname(model_probabilities)) %>%
    arrange(-`Posterior model probability`)
  
  model_probabilities <- model_probabilities %>%
    mutate(`% of top model probability` = round(100 * `Posterior model probability` / `Posterior model probability`[1], 1),
           `% of top model probability` = c("\\-", `% of top model probability`[2:length(`% of top model probability`)]))
  
 # saveRDS(model_list_expt2[[1]], "output/full_model_expt2.rds")
  saveRDS(model_list_expt2[[top_index]], "output/top_model_expt2.rds")
  saveRDS(model_probabilities, "output/model_probabilities_expt2.rds")
}

top_model <- readRDS("output/top_model_expt2.rds")
model_probabilities <- readRDS("output/model_probabilities_expt2.rds")

Posterior model probabilites

model_probabilities %>% 
  kable(digits = 3) %>% kable_styling()

Fixed effects from the top model

The model output suggests that the third possible outcome (i.e. being forced out of the hive; denoted as mu3 in the table) was more common in the ‘LPS CHCs’ treatment group. See Hypothesis testing and effect sizes below for more useful statistical results.

get_fixed_effects_with_p_values(top_model) %>% 
  kable(digits = 3) %>% kable_styling()

Derive prediction from the posterior

new <- expand.grid(treatment = levels(data_for_categorical_model$treatment), 
                   observation_time_minutes = 120)

preds <- fitted(top_model, newdata = new, summary = FALSE)
dimnames(preds) <- list(NULL, new[,1], NULL)

plotting_data <- rbind(
  as.data.frame(preds[,, 1]) %>% mutate(outcome = "Stayed inside the hive", posterior_sample = 1:n()),
  as.data.frame(preds[,, 2]) %>% mutate(outcome = "Left voluntarily", posterior_sample = 1:n()),
  as.data.frame(preds[,, 3]) %>% mutate(outcome = "Forced out", posterior_sample = 1:n())) %>%
  gather(treatment, prop, `Ringer CHC`, `LPS CHC`) %>%
  mutate(treatment = factor(treatment, c("Ringer CHC", "LPS CHC")),
         outcome = factor(outcome, c("Stayed inside the hive", "Left voluntarily", "Forced out"))) %>%
  as_tibble() %>% arrange(treatment, outcome) 

Make Figure 2

ggplot(plotting_data, aes(treatment, 100 * prop)) + 
  geom_eye(aes(fill = treatment)) + 
  scale_fill_brewer(palette = "Pastel1", direction = -1) + 
  facet_wrap( ~ outcome, scales = "free_y") + 
  ylab("% bees leaving the hive this way") + 
  xlab("Mode of exit") + 
  theme_bw() + 
  theme(legend.position = "none") 

% bees exiting the hive

my_summary <- function(df, columns, outcome) {
  lapply(columns, function(x){
    
    p <- df %>% pull(!! x) %>% bayestestR::p_direction() %>% as.numeric()
    p <- (100 - p) / 100
    
    df %>% pull(!! x) %>% posterior_summary() %>% as_tibble() %>% 
      mutate(p=p, Outcome = outcome, Metric = x) %>% select(Outcome, Metric, everything())
  }) %>% do.call("rbind", .)
}

stats_table <- rbind(
  plotting_data %>%
    filter(outcome == "Stayed inside the hive") %>%
    mutate(prop = 1 - prop) %>%
    spread(treatment, prop) %>%
    mutate(`Absolute difference in % bees exiting the hive` = 100 * (`LPS CHC` - `Ringer CHC`),
           `Log10 odds ratio` = log10(`LPS CHC` / `Ringer CHC`)) %>%
    my_summary(c("Absolute difference in % bees exiting the hive", 
                 "Log10 odds ratio"),
               outcome = "Stayed inside the hive") %>%
    mutate(p = c(" ", format(round(p[2], 4), nsmall = 4))),
  
  plotting_data %>%
    filter(outcome == "Left voluntarily") %>%
    spread(treatment, prop) %>%
    mutate(`Absolute difference in % bees leaving voluntarily` = 100 * (`LPS CHC` - `Ringer CHC`),
           `Log10 odds ratio` = log10(`LPS CHC` / `Ringer CHC`)) %>%
    my_summary(c("Absolute difference in % bees leaving voluntarily", 
                 "Log10 odds ratio"),
               outcome = "Left voluntarily") %>%
    mutate(p = c(" ", format(round(p[2], 4), nsmall = 4))),
  
  plotting_data %>%
    filter(outcome == "Forced out") %>%
    spread(treatment, prop) %>%
    mutate(`Absolute difference in % bees leaving voluntarily` = 100 * (`LPS CHC` - `Ringer CHC`),
           `Log10 odds ratio` = log10(`LPS CHC` / `Ringer CHC`)) %>%
    my_summary(c("Absolute difference in % bees leaving voluntarily", 
                 "Log10 odds ratio"),
               outcome = "Forced out") %>%
    mutate(p = c(" ", format(round(p[2], 4), nsmall = 4))) 
) 

stats_table[c(2,4,6), 1] <- " "

stats_table %>%
  mutate(` ` = ifelse(p < 0.05, "\\*", ""),
         ` ` = replace(` `, p < 0.01, "**"),
         ` ` = replace(` `, p < 0.001, "***"),
         ` ` = replace(` `, p == " ", "")) %>%
  kable(digits = 3) %>% kable_styling() %>% 
  row_spec(c(0,2,4,6), extra_css = "border-bottom: solid;")

Session information

sessionInfo()

R version 3.5.1 (2018-07-02)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS High Sierra 10.13.6

Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib

locale:
[1] en_AU.UTF-8/en_AU.UTF-8/en_AU.UTF-8/C/en_AU.UTF-8/en_AU.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] tidybayes_1.1.0  bayestestR_0.2.2 kableExtra_0.9.0 gridExtra_2.3   
 [5] forcats_0.4.0    stringr_1.4.0    dplyr_0.8.3      purrr_0.3.2     
 [9] readr_1.1.1      tidyr_0.8.2      tibble_2.1.3     ggplot2_3.1.0   
[13] tidyverse_1.2.1  bayesplot_1.6.0  brms_2.10.0      Rcpp_1.0.2      

loaded via a namespace (and not attached):
  [1] colorspace_1.3-2          ggridges_0.5.0           
  [3] rsconnect_0.8.8           rprojroot_1.3-2          
  [5] ggstance_0.3.1            markdown_1.0             
  [7] base64enc_0.1-3           fs_1.3.1                 
  [9] rstudioapi_0.10           rstan_2.19.2             
 [11] svUnit_0.7-12             DT_0.4                   
 [13] mvtnorm_1.0-11            lubridate_1.7.4          
 [15] xml2_1.2.0                bridgesampling_0.4-0     
 [17] knitr_1.23                shinythemes_1.1.1        
 [19] jsonlite_1.6              workflowr_1.4.0          
 [21] broom_0.5.0               shiny_1.3.2              
 [23] compiler_3.5.1            httr_1.4.0               
 [25] backports_1.1.2           assertthat_0.2.1         
 [27] Matrix_1.2-14             lazyeval_0.2.2           
 [29] cli_1.1.0                 later_0.8.0              
 [31] htmltools_0.3.6           prettyunits_1.0.2        
 [33] tools_3.5.1               igraph_1.2.1             
 [35] coda_0.19-2               gtable_0.2.0             
 [37] glue_1.3.1.9000           reshape2_1.4.3           
 [39] cellranger_1.1.0          nlme_3.1-137             
 [41] crosstalk_1.0.0           insight_0.3.0            
 [43] xfun_0.8                  ps_1.3.0                 
 [45] rvest_0.3.2               mime_0.7                 
 [47] miniUI_0.1.1.1            gtools_3.8.1             
 [49] zoo_1.8-3                 scales_1.0.0             
 [51] colourpicker_1.0          hms_0.4.2                
 [53] promises_1.0.1            Brobdingnag_1.2-5        
 [55] parallel_3.5.1            inline_0.3.15            
 [57] shinystan_2.5.0           RColorBrewer_1.1-2       
 [59] yaml_2.2.0                loo_2.1.0                
 [61] StanHeaders_2.19.0        stringi_1.4.3            
 [63] highr_0.8                 dygraphs_1.1.1.6         
 [65] pkgbuild_1.0.2            rlang_0.4.0              
 [67] pkgconfig_2.0.2           matrixStats_0.54.0       
 [69] evaluate_0.14             lattice_0.20-35          
 [71] rstantools_2.0.0          htmlwidgets_1.3          
 [73] labeling_0.3              processx_3.2.1           
 [75] tidyselect_0.2.5          plyr_1.8.4               
 [77] magrittr_1.5              R6_2.4.0                 
 [79] pillar_1.3.1.9000         haven_1.1.2              
 [81] whisker_0.3-2             withr_2.1.2              
 [83] xts_0.11-0                abind_1.4-5              
 [85] modelr_0.1.2              crayon_1.3.4             
 [87] arrayhelpers_1.0-20160527 rmarkdown_1.13           
 [89] grid_3.5.1                readxl_1.1.0             
 [91] callr_2.0.4               git2r_0.23.0             
 [93] threejs_0.3.1             digest_0.6.20            
 [95] xtable_1.8-4              httpuv_1.5.1             
 [97] stats4_3.5.1              munsell_0.5.0            
 [99] viridisLite_0.3.0         shinyjs_1.0