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", colour = "grey15") + 
  facet_wrap( ~ hive, nrow = 1, scales = "free_y") + 
  scale_fill_brewer(palette = "Set3", name = "Treatment") + 
  xlab("Outcome") + ylab("Number of bees") + 
  theme(axis.text.x = element_text(angle = 45, hjust=1),
        legend.position = "top")

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_expt1.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_expt1 <- 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_expt1[[1]], 
                                   model_list_expt1[[2]], 
                                   model_list_expt1[[3]], 
                                   model_list_expt1[[4]], 
                                   model_list_expt1[[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_expt1[[1]], "output/full_model_expt1.rds")
  saveRDS(model_list_expt1[[3]], "output/trt_model_expt1.rds")
  saveRDS(model_list_expt1[[top_index]], "output/top_model_expt1.rds")
  saveRDS(model_probabilities, "output/model_probabilities_expt1.rds")
}

top_model <- readRDS("output/top_model_expt1.rds")
full_model <- readRDS("output/full_model_expt1.rds")
model_probabilities <- readRDS("output/model_probabilities_expt1.rds")

Posterior model probabilites

This table shows the estimated probability that each of the 5 models we compared is the top model within this set of models. Note that the predictor variable treatment is not present in either of the top 2 models, but that the top models are not much more likely than the remaining models.

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

Fixed effects from the full model

The model output is not especially informative on its own, but is included here for completeness. See Hypothesis testing and effect sizes below for more useful statistical results.

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

Derive prediction from the posterior

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

preds <- fitted(full_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, `Intact control`, Ringers, `Heat-treated LPS`, LPS) %>% 
  mutate(outcome = factor(outcome, c("Stayed inside the hive", "Left voluntarily", "Forced out"))) %>%
  as_tibble() %>% arrange(treatment, outcome) 

stats_data <- plotting_data
plotting_data$treatment[plotting_data$treatment == "Intact control"] <- "Intact\ncontrol"
plotting_data$treatment[plotting_data$treatment == "Heat-treated LPS"] <- "Heat-treated\nLPS"
plotting_data$treatment <- factor(plotting_data$treatment, c("Intact\ncontrol", "Ringers", "Heat-treated\nLPS", "LPS"))

Make Figure 1

make_plot <- function(dat, lim){
  dat %>%
    ggplot(aes(treatment, 100 * prop)) + 
    geom_eye(aes(fill = treatment), size_range = c(0.5, 1.1), .width = c(0.5, 0.95)) + 
    scale_fill_brewer(palette = "Pastel2", direction = -1) + 
    facet_wrap( ~ outcome, scales = "free_y") + 
    coord_cartesian(ylim = lim) +
    ylab(NULL) + xlab(NULL) + 
    theme_bw() + 
    theme(legend.position = "none") 
}

grid.arrange(
  plotting_data %>% filter(outcome == "Stayed inside the hive") %>% make_plot(c(75, 100)),
  plotting_data %>% filter(outcome == "Left voluntarily") %>% make_plot(c(0, 10)),
  plotting_data %>% filter(outcome == "Forced out") %>% make_plot(c(0, 25)),
  left = "% bees leaving the hive this way", bottom = "Treatment", nrow = 1)

Figure 1: This figure shows the estimated % bees from each treatment group that either stayed inside the hive (left), left the hive voluntarily (middle), or were forced out by their nestmates (right). The points show the median of the posterior distribution of estimates from the model, while the thick and thin error bars show the 50% and 95% credible intervals of this posterior distribution, respectively. The coloured area shows the whole posterior. Bees in the ‘intact control’ group were more likely to stay inside the hive, and were less likely to be forced out, than bees in either heat-treated LPS or LPS treatments (see Table XX).

% bees exiting the hive

my_summary <- function(df, columns) {
  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) %>% mutate(Metric = x) %>% select(Metric, everything())
  }) %>% do.call("rbind", .)
}


make_stats_table <- function(dat, groupA, groupB, 
                             comparison, 
                             metric = "Absolute difference in % bees exiting the hive"){
  output <- dat %>%
    spread(treatment, prop) %>%
    mutate(
       metric_here = 100 * (!! enquo(groupA) - !! enquo(groupB)), 
      `Log10 odds ratio` = log10(!! enquo(groupA) / !! enquo(groupB))) %>%
    my_summary(c("metric_here", "Log10 odds ratio")) %>%
    mutate(p = c(" ", format(round(p[2], 4), nsmall = 4)),
           Comparison = comparison) %>%
    select(Comparison, everything()) %>%
    mutate(Metric = replace(Metric, Metric == "metric_here", metric))
  
  names(output)[names(output) == "metric_here"] <- metric
  output 
}

stats_table <- rbind(
  stats_data %>%
    filter(outcome == "Stayed inside the hive") %>%
    make_stats_table(`Intact control`, `Heat-treated LPS`, "Intact control vs Heat-treated LPS"),
  
  stats_data %>%
    filter(outcome == "Stayed inside the hive") %>%
    make_stats_table(`Ringers`, `Heat-treated LPS`, "Ringers control vs Heat-treated LPS"),
  
  stats_data %>%
    filter(outcome == "Stayed inside the hive") %>%
    make_stats_table(`Intact control`, `LPS`, "Intact control vs LPS"),
  
  stats_data %>%
    filter(outcome == "Stayed inside the hive") %>%
    make_stats_table(`Ringers`, `LPS`, "Ringers control vs LPS")
) %>% as_tibble()

stats_table[c(2,4,6,8), 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,8), extra_css = "border-bottom: solid;")

% bees that were forced out

stats_table2 <- rbind(
  stats_data %>%
    filter(outcome == "Forced out") %>%
    make_stats_table(`Intact control`, `Heat-treated LPS`, 
                     "Intact control vs Heat-treated LPS", 
                     metric = "Absolute difference in % bees that were forced out"),
  
  stats_data %>%
    filter(outcome == "Forced out") %>%
    make_stats_table(`Ringers`, `Heat-treated LPS`, 
                     "Ringers control vs Heat-treated LPS", 
                     metric = "Absolute difference in % bees that were forced out"),
  
  stats_data %>%
    filter(outcome == "Forced out") %>%
    make_stats_table(`Intact control`, `LPS`, 
                     "Intact control vs LPS", 
                     metric = "Absolute difference in % bees that were forced out"),
  
  stats_data %>%
    filter(outcome == "Forced out") %>%
    make_stats_table(`Ringers`, `LPS`, 
                     "Ringers control vs LPS", 
                     metric = "Absolute difference in % bees that were forced out")
) %>% as_tibble()

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

stats_table2 %>%
  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,8), extra_css = "border-bottom: solid;")

% bees that left voluntarily

stats_table3 <- rbind(
  stats_data %>%
    filter(outcome == "Left voluntarily") %>%
    make_stats_table(`Intact control`, `Heat-treated LPS`, 
                     "Intact control vs Heat-treated LPS", 
                     metric = "Absolute difference in % bees that left voluntarily"),
  
  stats_data %>%
    filter(outcome == "Left voluntarily") %>%
    make_stats_table(`Ringers`, `Heat-treated LPS`, 
                     "Ringers control vs Heat-treated LPS", 
                     metric = "Absolute difference in % bees that left voluntarily"),
  
  stats_data %>%
    filter(outcome == "Left voluntarily") %>%
    make_stats_table(`Intact control`, `LPS`, 
                     "Intact control vs LPS", 
                     metric = "Absolute difference in % bees that left voluntarily"),
  
  stats_data %>%
    filter(outcome == "Left voluntarily") %>%
    make_stats_table(`Ringers`, `LPS`, 
                     "Ringers control vs LPS", 
                     metric = "Absolute difference in % bees that left voluntarily")
) %>% as_tibble()

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

stats_table3 %>%
  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,8), 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