statistics

Load R libraries

library(tidyverse)
library(brms)
library(bayestestR)
library(kableExtra)
library(ggbeeswarm)
library(RColorBrewer)
library(showtext)
font_add_google(name = "Lato", family = "Lato", regular.wt = 400, bold.wt = 700)
showtext_auto()
options(stringsAsFactors = FALSE)

SE <- function(x) sd(x) / sqrt(length(x))

Load the data

fitness_data <- read_csv("data/SR_fitness_data.csv") %>% 
  filter(!is.na(genotype)) %>%
  rename(body_size = `Body size`,
         female_age = `F age`) %>%
  mutate(genotype = factor(genotype, levels = c("STST", "SRST", "SRSR")))

Warning: Missing column names filled in: 'X12' [12]

Parsed with column specification:
cols(
  ID = col_integer(),
  Isoline = col_character(),
  genotype = col_character(),
  `genotype pooled 2` = col_integer(),
  `genotype pooled 1` = col_integer(),
  `genotype pooled` = col_integer(),
  Gen = col_integer(),
  Block = col_integer(),
  `F age` = col_integer(),
  `Body size` = col_double(),
  offspring = col_integer(),
  X12 = col_character()
)

Make a table of summary statistics and sample sizes

Here, we calculate the mean offspring produced by females from each of the three genotypes (STST, SRST, and SRSR), either within each isoline or across all the isolines. We also calculate the % females that failed to produce any offspring, and provide sample size information.

means_by_isoline <- fitness_data %>%
  group_by(genotype, Isoline) %>%
  summarise(
    Number_of_females_measured = n(),
    Mean_offspring_per_female = mean(offspring),
    SE = SE(offspring),
    n_females_producing_offspring = sum(offspring != 0),
    Percent_females_producing_offspring = 100 * n_females_producing_offspring / n()) 

means <- fitness_data %>%
  mutate(Isoline = "Across all isolines") %>%
  group_by(genotype, Isoline) %>%
  summarise(
    Number_of_females_measured = n(),
    Mean_offspring_per_female = mean(offspring),
    SE = SE(offspring),
    n_females_producing_offspring = sum(offspring != 0),
    Percent_females_producing_offspring = 100 * n_females_producing_offspring / n()) 


bind_rows(means_by_isoline, means) %>%
  rename_all(function(x) gsub("_", " ", x)) %>%
  rename_all(function(x) gsub("Percent", "%", x)) %>%
  rename(Genotype = genotype) %>%
  kable(digits = 2) %>% kable_styling()

Run the Bayesian hurdle model

The model assumes that the response variable, offspring number, is the result of a ‘hurdle’ process. Essentially this means that the model consists of two sub-models: one controlling the probability that offspring number is non-zero, and one controlling the number of offspring produced provided that more than zero are produced (we assume that offspring number follows a negative binomial distribution, because this improved model fit relative to the simpler hurdle-Poisson model).

We assume that the parameters controlling both the hurdle and the distribution of non-zero values are affected by four fixed effects (the female’s genotype: STST, SRST, or SRSR), her isoline, the female’s age, and the interaction between genotype and isoline. We also fit two random effects: isoline, and experimental block. All fixed effects were assumed to have a prior distribution following a normal distribution with mean 0 and SD = 5.

if(!file.exists("output/brms_model.rds")){

  model_formula <- bf(
    offspring ~ genotype * Isoline + female_age + (1 | Block), # formula for the mean
    hu        ~ genotype * Isoline + female_age + (1 | Block)  # formula for the hurdle
  )
  
  model_prior <- c(set_prior("normal(0, 5)", class = "b"),
                   set_prior("normal(0, 5)", class = "b", dpar = "hu"))
  
  # imputed_data <- fitness_data %>%
  #   select(offspring, genotype, body_size, female_age, Isoline, Block) %>%
  #   mice(m = 5, print = FALSE) 
  
  model <- brm(model_formula,
                        family = "hurdle_negbinomial",
                        chains = 4, cores = 1, iter = 4000, inits = 0, seed = 12345,
                        control = list(adapt_delta = 0.999, max_treedepth = 15),
                        prior = model_prior,
                        data = fitness_data)
  saveRDS(model, file = "output/brms_model.rds")
} else model <- readRDS("output/brms_model.rds")

Perform a posterior predictive check

The idea behind posterior predictive checking is that if our model is a good fit, then we should be able to use it to generate data that looks a lot like the data we observed. Here, we see 10 draws from the ‘posterior predictive distribution’ (thin lines), which indeed look quite similar to the distribution of the real data (thick line), suggested that our model approximates the processes that generated the real data well enough for reliable inference.

pp_check(model, type = "dens_overlay")

Using 10 posterior samples for ppc type 'dens_overlay' by default.

Inspect the model’s parameter estimates

bayesian_p_values <- as.data.frame(p_direction(model)) %>% 
      mutate(pd = (100 - pd) / 100,
             Parameter = gsub("[.]", ":", gsub("b_", "", Parameter)))

random <- summary(model)$random # Finish

summary(model)$fixed %>% as.data.frame() %>% 
  rownames_to_column("Parameter") %>%
  left_join(bayesian_p_values, by = "Parameter") %>%
  rename(p = pd) %>% arrange(grepl("hu_", Parameter)) %>%
  mutate(Parameter = gsub("hu_", "Hurdle - ", Parameter),
         Estimate =  format(round(Estimate, 3), nsmall = 3),
         Est.Error =  format(round(Est.Error, 3), nsmall = 3),
         ` ` = ifelse(p < 0.05, "*", ""),
         p = format(round(p, 4), nsmall = 4),
         Rhat = format(round(Rhat, 3), nsmall = 3),
         `l-95% CI` = format(round(`l-95% CI`, 3), nsmall = 3),
         `u-95% CI` = format(round(`u-95% CI`, 3), nsmall = 3),
         Eff.Sample = round(Eff.Sample, 0)
         ) 

                            Parameter Estimate Est.Error l-95% CI u-95% CI
1                           Intercept    4.049     0.269    3.520    4.585
2                        genotypeSRST    0.366     0.169    0.029    0.698
3                        genotypeSRSR   -0.423     0.186   -0.790   -0.065
4                        IsolineLew17   -0.061     0.164   -0.378    0.261
5                        IsolineSloB3    0.287     0.167   -0.048    0.612
6                        IsolineSloB7    0.163     0.171   -0.170    0.496
7                          female_age    0.023     0.058   -0.091    0.138
8           genotypeSRST:IsolineLew17   -0.225     0.242   -0.690    0.261
9           genotypeSRSR:IsolineLew17   -0.017     0.255   -0.513    0.488
10          genotypeSRST:IsolineSloB3   -0.687     0.253   -1.187   -0.191
11          genotypeSRSR:IsolineSloB3   -0.849     0.268   -1.361   -0.314
12          genotypeSRST:IsolineSloB7   -0.547     0.240   -1.015   -0.069
13          genotypeSRSR:IsolineSloB7   -0.340     0.263   -0.854    0.176
14                 Hurdle - Intercept   -1.498     1.190   -3.834    0.816
15              Hurdle - genotypeSRST    1.473     0.799    0.003    3.142
16              Hurdle - genotypeSRSR    2.186     0.789    0.750    3.854
17              Hurdle - IsolineLew17   -1.045     1.206   -3.607    1.140
18              Hurdle - IsolineSloB3   -0.925     1.195   -3.408    1.287
19              Hurdle - IsolineSloB7   -0.666     1.197   -3.150    1.504
20                Hurdle - female_age   -0.327     0.237   -0.796    0.131
21 Hurdle - genotypeSRST:IsolineLew17    0.564     1.344   -1.982    3.384
22 Hurdle - genotypeSRSR:IsolineLew17   -0.833     1.399   -3.483    2.078
23 Hurdle - genotypeSRST:IsolineSloB3    0.851     1.343   -1.656    3.658
24 Hurdle - genotypeSRSR:IsolineSloB3    1.078     1.307   -1.358    3.753
25 Hurdle - genotypeSRST:IsolineSloB7   -0.590     1.402   -3.270    2.293
26 Hurdle - genotypeSRSR:IsolineSloB7    0.186     1.291   -2.219    2.827
   Eff.Sample  Rhat      p  
1        6153 1.000 0.0000 *
2        3221 1.001 0.0166 *
3        3580 1.001 0.0111 *
4        4053 1.000 0.3622  
5        3999 1.000 0.0424 *
6        3623 1.000 0.1717  
7        9705 1.000 0.3435  
8        3965 1.000 0.1704  
9        4135 1.001 0.4712  
10       4235 1.000 0.0045 *
11       4732 1.001 0.0009 *
12       3817 1.000 0.0115 *
13       4039 1.001 0.0988  
14       5225 1.000 0.1005  
15       3454 1.000 0.0250 *
16       3472 1.000 0.0009 *
17       3525 1.000 0.1933  
18       2814 1.000 0.2150  
19       3975 1.000 0.3049  
20      10966 1.000 0.0831  
21       3611 1.000 0.3439  
22       3812 1.000 0.2704  
23       3169 1.000 0.2689  
24       2985 1.000 0.2051  
25       4095 1.000 0.3333  
26       3847 1.000 0.4508  

Plot the estimated means for each genotype

new <- fitness_data %>% select(genotype, Isoline, body_size, female_age) %>%
  mutate(body_size  = mean(body_size, na.rm = TRUE),
         female_age = mean(female_age)) %>% 
  distinct()

predicted_mean <- data.frame(new, fitted(model, newdata = new, re_formula = NA)) %>% 
  mutate(facet = "A. Mean offspring production")
predicted_mean_when_fertile <- data.frame(new, fitted(model, newdata = new, dpar = "mu", re_formula = NA)) %>%
  mutate(facet = "B. Mean offspring production\n(excluding infertile females)")
predicted_prop_infertile <- data.frame(new, fitted(model, newdata = new, dpar = "hu", re_formula = NA)) %>% 
  mutate(facet = "C. % fertile females",
         Estimate = 100 * (1 - Estimate), Q2.5 = (1 - Q2.5) * 100, Q97.5 = (1 - Q97.5) * 100)


posterior_means <- 
  data.frame(new, as.data.frame(t(fitted(model, newdata = new, re_formula = NA, summary = FALSE)))) %>%
  select(-body_size, -female_age, -Isoline) %>%
  group_by(genotype) %>%
  summarise_all(mean) %>% select(-genotype) %>% t() %>% as.data.frame() %>%
  rename(STST = V1, SRST = V2, SRSR = V3)
  
predicted_mean_all_iso <- data.frame(
  genotype = c("STST", "SRST", "SRSR"), Isoline = "All", 
  posterior_means %>% lapply(posterior_summary) %>% do.call("rbind", .),
  facet = "A. Mean offspring production") 

posterior_means_when_fertile <- 
  data.frame(new, as.data.frame(t(fitted(model, newdata = new, dpar = "mu", re_formula = NA, summary = FALSE)))) %>%
  select(-body_size, -female_age, -Isoline) %>%
  group_by(genotype) %>%
  summarise_all(mean) %>% select(-genotype) %>% t() %>% as.data.frame() %>%
  rename(STST = V1, SRST = V2, SRSR = V3)
  
predicted_mean_when_fertile_all_iso <- data.frame(
  genotype = c("STST", "SRST", "SRSR"), Isoline = "All", 
  posterior_means_when_fertile %>% lapply(posterior_summary) %>% do.call("rbind", .),
  facet = "B. Mean offspring production\n(excluding infertile females)") 

posterior_means_prop_infertile <- 
  data.frame(new, as.data.frame(t(fitted(model, newdata = new, dpar = "hu", re_formula = NA, summary = FALSE)))) %>%
  select(-body_size, -female_age, -Isoline) %>%
  group_by(genotype) %>%
  summarise_all(~ 100 * mean(.x)) %>% select(-genotype) %>% t() %>% as.data.frame() %>%
  rename(STST = V1, SRST = V2, SRSR = V3)
  
predicted_prop_infertile_all_iso <- data.frame(
  genotype = c("STST", "SRST", "SRSR"), Isoline = "All", 
  posterior_means_prop_infertile %>% lapply(posterior_summary) %>% do.call("rbind", .),
  facet = "C. % fertile females") 

preds_figure1 <- bind_rows(predicted_mean_all_iso,
                   predicted_mean_when_fertile_all_iso,
                   predicted_prop_infertile_all_iso) %>%
  mutate(genotype = factor(genotype, levels = c("STST", "SRST", "SRSR")))

preds_figure2 <- bind_rows(predicted_mean,
                         predicted_mean_when_fertile,
                         predicted_prop_infertile) %>%
  mutate(facet = factor(facet, levels = unique(facet)),
         genotype = factor(genotype, levels = c("STST", "SRST", "SRSR"))) %>%
  select(-body_size, -female_age)


beeswarm_points <- bind_rows(
  fitness_data %>% mutate(facet = "A. Mean offspring production"),
  fitness_data %>% filter(offspring != 0) %>% mutate(facet = "B. Mean offspring production\n(excluding infertile females)")) %>% 
  mutate(Fertility = ifelse(offspring == 0, "Sterile", "Fertile"),
         genotype  = factor(genotype, levels = c("STST", "SRST", "SRSR"))) %>%
    rename(Estimate = offspring) 
  

pos1 <- position_nudge(x = -0.17)
pal <- c(brewer.pal(4, "RdPu")[2], brewer.pal(7, "Purples")[5])

figure_1 <- preds_figure1 %>%
  ggplot(aes(genotype, Estimate)) + 
  geom_errorbar(aes(ymin = Q2.5, ymax = Q97.5), colour = "grey20", position = pos1, size = .8, width = 0.1) + 
  geom_beeswarm(data = beeswarm_points, aes(colour = Fertility),
                size = .7, alpha = 0.6) + 
  geom_point(size = 3.1, pch = 21, colour = "black", position = pos1, fill = "grey20") + 
  scale_colour_manual(values = pal) + 
  facet_wrap(~facet, scale = "free_y") + 
  labs(y = "Posterior estimate \u00B1 95% CIs", x = "Genotype") + 
  theme_bw() + 
  theme(strip.background = element_blank(),
        text = element_text(family = "Lato", size = 12),
        panel.grid.major.x = element_blank(), 
        strip.text = element_text(hjust = 0))


dodge <- position_dodge(0.66)
figure_2 <- preds_figure2 %>%
  ggplot(aes(genotype, Estimate, fill = Isoline)) + 
  geom_errorbar(aes(ymin = Q2.5, ymax = Q97.5), size = .7, width = 0.3, colour = "grey40", position = dodge) + 
  geom_point(size = 3.1, pch = 21, colour = "black", position = dodge) + 
  facet_wrap(~facet, scale = "free_y") + 
  scale_fill_brewer(palette = "Pastel1") +
  labs(y = "Posterior estimate \u00B1 95% CIs", x = "Genotype") + 
  theme_bw() + 
  theme(strip.background = element_blank(),
        text = element_text(family = "Lato", size = 12),
        panel.grid.major.x = element_blank(), 
        strip.text = element_text(hjust = 0))

figure_1 %>% ggsave(filename = "figures/figure_1.pdf", width = 9, height = 4)
figure_2 %>% ggsave(filename = "figures/figure_2.pdf", width = 9, height = 4)

figure_1

figure_2

Calculated pairwise differences between genotypes

compare_means <- function(mean1, mean2, posterior){
  difference <- posterior[, mean2] - posterior[, mean1]
  relative <- median(posterior[, mean2] / posterior[, mean1])
  p_value <- as.numeric(p_direction(difference)) / 100
  as_tibble(posterior_summary(as.mcmc(difference))) %>%
    mutate(Comparison = paste(mean1, mean2, sep = " \u2192 "),
           `% difference` = format(round(100 * relative, 1), nsmall = 1),
           `95% CIs` = paste(" (", format(round(Q2.5, 1), nsmall = 1), " to ", format(round(Q97.5, 1), nsmall = 1), "%)", sep = ""),
           `Difference in means` = paste(format(round(Estimate, 2), nsmall = 2), `95% CIs`, sep = ""),
           `Fitness trait` = NA,
           p = p_value) %>%
    select( -Q2.5, -Q97.5) %>%
    select(Comparison, `Fitness trait`, `Difference in means`, Est.Error, `% difference`, p)
}

table_of_contrasts <- bind_rows(
  compare_means("STST", "SRST", posterior_means),
  compare_means("STST", "SRSR", posterior_means),
  compare_means("SRST", "SRSR", posterior_means),
  compare_means("STST", "SRST", posterior_means_when_fertile),
  compare_means("STST", "SRSR", posterior_means_when_fertile),
  compare_means("SRST", "SRSR", posterior_means_when_fertile),
  compare_means("STST", "SRST", posterior_means_prop_infertile),
  compare_means("STST", "SRSR", posterior_means_prop_infertile),
  compare_means("SRST", "SRSR", posterior_means_prop_infertile)
) %>% mutate(`Fitness trait` = rep(c("Mean offspring production",
                                     "Mean offspring production (excluding infertile females)",
                                     "% fertile females"), each = 3)) %>%
  mutate(Est.Error = format(round(Est.Error, 2), nsmall = 2),
         ` ` = ifelse(p < 0.05, "*", " "),
         p = format(round(p, 4), nsmall = 4))

table_of_contrasts %>%
  kable() %>% kable_styling()

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] showtext_0.5-1     showtextdb_2.0     sysfonts_0.7.2    
 [4] RColorBrewer_1.1-2 ggbeeswarm_0.6.0   kableExtra_0.9.0  
 [7] bayestestR_0.2.2   brms_2.8.0         Rcpp_1.0.1        
[10] forcats_0.4.0      stringr_1.4.0      dplyr_0.8.0.1     
[13] purrr_0.3.2        readr_1.1.1        tidyr_0.8.2       
[16] tibble_2.0.99.9000 ggplot2_3.1.0      tidyverse_1.2.1   

loaded via a namespace (and not attached):
  [1] colorspace_1.3-2     ggridges_0.5.0       rsconnect_0.8.8     
  [4] rprojroot_1.3-2      markdown_0.9         base64enc_0.1-3     
  [7] fs_1.3.1             rstudioapi_0.10      rstan_2.18.2        
 [10] DT_0.4               mvtnorm_1.0-8        lubridate_1.7.4     
 [13] xml2_1.2.0           bridgesampling_0.4-0 knitr_1.22          
 [16] shinythemes_1.1.1    bayesplot_1.6.0      jsonlite_1.6        
 [19] workflowr_1.3.0      broom_0.5.0          shiny_1.3.2         
 [22] compiler_3.5.1       httr_1.4.0           backports_1.1.2     
 [25] assertthat_0.2.1     Matrix_1.2-14        lazyeval_0.2.2      
 [28] cli_1.1.0            later_0.8.0          htmltools_0.3.6     
 [31] prettyunits_1.0.2    tools_3.5.1          igraph_1.2.1        
 [34] coda_0.19-2          gtable_0.2.0         glue_1.3.1.9000     
 [37] reshape2_1.4.3       cellranger_1.1.0     nlme_3.1-137        
 [40] crosstalk_1.0.0      insight_0.3.0        xfun_0.6            
 [43] ps_1.3.0             rvest_0.3.2          mime_0.6            
 [46] miniUI_0.1.1.1       gtools_3.8.1         zoo_1.8-3           
 [49] scales_1.0.0         colourpicker_1.0     hms_0.4.2           
 [52] promises_1.0.1       Brobdingnag_1.2-5    parallel_3.5.1      
 [55] inline_0.3.15        shinystan_2.5.0      curl_3.3            
 [58] yaml_2.2.0           gridExtra_2.3        loo_2.1.0           
 [61] StanHeaders_2.18.0   stringi_1.4.3        highr_0.8           
 [64] dygraphs_1.1.1.6     pkgbuild_1.0.2       rlang_0.3.4         
 [67] pkgconfig_2.0.2      matrixStats_0.54.0   evaluate_0.13       
 [70] lattice_0.20-35      labeling_0.3         rstantools_1.5.0    
 [73] htmlwidgets_1.2      tidyselect_0.2.5     processx_3.2.1      
 [76] plyr_1.8.4           magrittr_1.5         R6_2.4.0            
 [79] pillar_1.3.1.9000    haven_1.1.2          withr_2.1.2         
 [82] xts_0.11-0           abind_1.4-5          modelr_0.1.2        
 [85] crayon_1.3.4         rmarkdown_1.10       grid_3.5.1          
 [88] readxl_1.1.0         callr_2.0.4          git2r_0.23.0        
 [91] threejs_0.3.1        digest_0.6.18        xtable_1.8-4        
 [94] httpuv_1.5.1         stats4_3.5.1         munsell_0.5.0       
 [97] beeswarm_0.2.3       viridisLite_0.3.0    vipor_0.4.5         
[100] shinyjs_1.0