Statistical analysis

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")){
 
  # The hurdle and the mean have the same set of predictors
  model_formula <- bf(
    offspring ~ genotype * Isoline + female_age + (1 | Block), 
    hu        ~ genotype * Isoline + female_age + (1 | Block)  
  )
  
  model_formula2 <- bf(
    offspring ~ genotype + Isoline + female_age + (1 | Block), 
    hu        ~ genotype + Isoline + female_age + (1 | Block)  
  )
    
  model_formula3 <- bf(
    offspring ~ genotype + female_age + (1 | Block), 
    hu        ~ genotype + female_age + (1 | Block)  
  )
  
  # This model added in response to a reviewer request. 
  # I don't think this model is much use here, since body size is a mediator variable;
  # see the diagram at the top of this page: https://en.wikipedia.org/wiki/Mediation_(statistics)
  # Moreover, many females have no body size measurement, so the sample size is smaller
  body_size_model <- bf(
    offspring ~ genotype + female_age + body_size + (1 | Block), 
    hu        ~ genotype + female_age + body_size + (1 | Block)  
  )
  
  
  # Find R^2 for a brms model, and its 95% CIs, and present neatly
  neat_R2 <- function(model){
    R2 <- bayes_R2(model) %>% round(2)
    paste(R2[1,1], " (95% CIs = ", R2[1,3], "-", R2[1,4], ")", sep = "")
  }
  
  # We set conservative, "regularising" priors - see McElreath's "Statistical Rethinking" textbook
  model_prior <- c(set_prior("normal(0, 3)", class = "b"),
                   set_prior("normal(0, 3)", class = "b", dpar = "hu"))
  
  full_model <- brm(model_formula,
                    family = "hurdle_negbinomial",
                    chains = 4, cores = 1, iter = 40000, inits = 0, seed = 12345,
                    control = list(adapt_delta = 0.9999, max_treedepth = 15),
                    save_all_pars = TRUE, 
                    prior = model_prior, 
                    data = fitness_data)
  
  no_interaction <- brm(model_formula2,
                        family = "hurdle_negbinomial",
                        chains = 4, cores = 1, iter = 40000, inits = 0, seed = 12345,
                        control = list(adapt_delta = 0.9999, max_treedepth = 15),
                        save_all_pars = TRUE, 
                        prior = model_prior, 
                        data = fitness_data)
  
  genotype_only_model <- brm(model_formula3,
                             family = "hurdle_negbinomial",
                             chains = 4, cores = 1, iter = 40000, inits = 0, seed = 12345,
                             control = list(adapt_delta = 0.9999, max_treedepth = 15),
                             save_all_pars = TRUE, 
                             prior = model_prior, 
                             data = fitness_data)
  
  # Added after peer review:
  body_size_model <- brm(body_size_model,
                             family = "hurdle_negbinomial",
                             chains = 4, cores = 1, iter = 40000, inits = 0, seed = 12345,
                             control = list(adapt_delta = 0.9999, max_treedepth = 15),
                             save_all_pars = TRUE, 
                             prior = model_prior, 
                             data = fitness_data)
  
  saveRDS(post_prob(full_model, no_interaction, genotype_only_model), 
          file = "output/model_comparison.rds")
  saveRDS(full_model, file = "output/full_model.rds")
  saveRDS(genotype_only_model, file = "output/genotype_only_model.rds")
  saveRDS(body_size_model, file = "output/body_size_model.rds")
  saveRDS(neat_R2(full_model), file = "output/R2_of_full_model.rds")
  saveRDS(neat_R2(genotype_only_model), file = "output/R2_of_genotype_only_model.rds")
  saveRDS(neat_R2(body_size_model), file = "output/R2_of_body_size_model.rds")
} else{
  full_model <- readRDS("output/full_model.rds")
  genotype_only_model <- readRDS("output/genotype_only_model.rds")
  model_probabilities <- readRDS("output/model_comparison.rds")
  body_size_model <- readRDS("output/body_size_model.rds")
}

Graphically verify the fit of the model using 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 a dataset which looks a lot like the dataset we actually observed. Here, we see 11 draws from the ‘posterior predictive distribution’ (pale blue), which indeed look quite similar to the distribution of the real data (dark blue), suggesting that our model is a good enough approximation of the true data-generating process for reliable inference.

pp_check(genotype_only_model, type = "hist", nsamples = 11, binwidth = 5)

Past versions of unnamed-chunk-4-1.png

Generate posterior predictions of the group means

Here, we estimate the mean for three measures of female fitness using the model, for each genotype (across all isolines) and for each genotype-isoline combination. The model adjusts for variation due to experimental block and female age.

make_figure_data <- function(by_isoline = FALSE){
  if(by_isoline){
    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()
    model <- full_model
    col_names <- paste(new$genotype, new$Isoline, sep = "~")
  } else {
    new <- fitness_data %>% 
      select(genotype, body_size, female_age) %>%
      mutate(body_size  = mean(body_size, na.rm = TRUE),
             female_age = mean(female_age)) %>% 
      distinct()
    model <- genotype_only_model
    col_names <- new$genotype
  }
  
  # Summarise the posterior (dots and CIs in Figure 1 or S1)
  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_fertile <- data.frame(new, fitted(model, newdata = new, dpar = "hu", re_formula = NA)) %>% 
    mutate(facet = "C. % fertile females",
           Estimate = 100 * (1 - Estimate), # Convert to percentage of fertile females, instead of *proportion* that are *in*fertile
           Q2.5 = (1 - Q2.5) * 100, 
           Q97.5 = (1 - Q97.5) * 100)
  
  summary_df <- bind_rows(predicted_mean,
                          predicted_mean_when_fertile,
                          predicted_prop_fertile) %>%
    mutate(genotype = factor(genotype, levels = c("STST", "SRST", "SRSR")))
  if(!by_isoline) summary_df <- summary_df %>% mutate(Isoline = "All isolines")
  
  # Posterior for facet A (overal progeny)
  posterior_means <- fitted(model, newdata = new, re_formula = NA, summary = FALSE) %>% as.data.frame()
  names(posterior_means) <- col_names
  posterior_facetA <- gather(posterior_means) %>% 
    mutate(facet = "A. Mean offspring production")
  
  # Posterior for facet B (excluding infertile females)
  posterior_means <- fitted(model, newdata = new, dpar = "mu", re_formula = NA, summary = FALSE) %>% as.data.frame()
  names(posterior_means) <- col_names
  posterior_facetB <- gather(posterior_means) %>% 
    mutate(facet = "B. Mean offspring production\n(excluding infertile females)")
  
  # Posterior for facet C (% infertile females)
  posterior_means <- fitted(model, newdata = new, dpar = "hu", re_formula = NA, summary = FALSE) %>% as.data.frame()
  names(posterior_means) <- col_names
  posterior_facetC <- gather(posterior_means) %>% 
    mutate(facet = "C. % fertile females")
  
  posterior_df <- bind_rows(
    posterior_facetA, posterior_facetB, posterior_facetC
  )
  
  if(by_isoline){
    posterior_df <- posterior_df %>%
      mutate(split = strsplit(key, split = "~"), 
             genotype = map_chr(split, ~ .x[1]),
             Isoline = map_chr(split, ~ .x[2])) %>% select(-key)
  } else {
    posterior_df <- posterior_df %>%
      rename(genotype = key)
  }
  
  posterior_df <- posterior_df %>% mutate(genotype = factor(genotype, levels = c("STST", "SRST", "SRSR")))
  
  list(summary_df, posterior_df)
}

figure1_data <- make_figure_data()
figureS1_data <- make_figure_data(by_isoline = TRUE)

Plot the posterior predictions of the group means

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) 
  
pal <- c("#6ca0dc", "#e34132")

figure_1 <- figure1_data[[1]] %>%
  ggplot(aes(genotype, Estimate)) + 
  geom_quasirandom(data = beeswarm_points, aes(colour = Fertility),
                size = .7, alpha = 0.6) + 
  geom_errorbar(aes(ymin = Q2.5, ymax = Q97.5), colour = "grey20", size = .8, width = 0) + 
  geom_point(size = 3.1, pch = 21, colour = "black", 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))

figure_S1 <- figureS1_data[[1]] %>%
  ggplot(aes(genotype, Estimate, fill = Isoline)) + 
  geom_errorbar(aes(ymin = Q2.5, ymax = Q97.5), size = .7, width = 0, colour = "grey40", position = position_dodge(0.7)) + 
  geom_point(size = 3.1, pch = 21, colour = "black", position = position_dodge(0.7)) + 
  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_S1 %>% ggsave(filename = "figures/figure_S1.pdf", width = 9, height = 4)

figure_1

Past versions of unnamed-chunk-10-1.png

Figure 1: The black points and error bars show the posterior estimates of the genotype means for A) offspring production, B) offspring production among the set of females that produced at least one offspring, and C) the percentage of females that produced offspring. The estimates are all derived from a single hurdle model which adjusts for variation due to female age and experimental block, and each estimate is the average across the four isolines (see Figure S1 for estimates split by isoline). The points show the raw values of offspring production for individual females, and are coloured purple for females that produced no offspring. The error bars show the 95% credible intervals on each estimate.

figure_S1

Figure S1: The same information as in Figure 1, except split by isoline.

Associated statistics for Figure 2

body_size_model_data <- fitness_data %>%
  filter(!is.na(body_size)) %>%
  mutate(body_size_scaled = as.numeric(scale(body_size)))

body_size_model_data %>%
  group_by(genotype) %>%
  summarise(`Mean wing vein length (mm)` = mean(body_size),
            SE = sd(body_size) / sqrt(n())) %>%
  kable(digits=3) %>% kable_styling()

body_size_model <- lmer(body_size_scaled ~ genotype + Isoline + (1 | Block), 
                      data = body_size_model_data)
summary(body_size_model)

Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: body_size_scaled ~ genotype + Isoline + (1 | Block)
   Data: body_size_model_data

REML criterion at convergence: 799.9

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.8321 -0.5805  0.1347  0.6665  3.6208 

Random effects:
 Groups   Name        Variance Std.Dev.
 Block    (Intercept) 0.0000   0.0000  
 Residual             0.6028   0.7764  
Number of obs: 338, groups:  Block, 6

Fixed effects:
               Estimate Std. Error        df t value Pr(>|t|)    
(Intercept)    -0.03817    0.09786 332.00000  -0.390   0.6967    
genotypeSRST    1.10351    0.10863 332.00000  10.158  < 2e-16 ***
genotypeSRSR    0.45672    0.10970 332.00000   4.163  4.0e-05 ***
IsolineLew 17  -1.10869    0.10595 332.00000 -10.464  < 2e-16 ***
IsolineSlo B3  -0.52049    0.12968 332.00000  -4.014  7.4e-05 ***
IsolineSlo B7  -0.27158    0.12367 332.00000  -2.196   0.0288 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) gnSRST gnSRSR IslL17 IslSB3
genotypSRST -0.573                            
genotypSRSR -0.574  0.556                     
IsolineLw17 -0.550  0.021  0.025              
IsolineSlB3 -0.272 -0.219 -0.226  0.396       
IsolineSlB7 -0.554  0.126  0.138  0.426  0.307

anova(body_size_model)

Type III Analysis of Variance Table with Satterthwaite's method
         Sum Sq Mean Sq NumDF DenDF F value    Pr(>F)    
genotype 64.133  32.067     2   332  53.194 < 2.2e-16 ***
Isoline  69.825  23.275     3   332  38.610 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

sample_sizes <- body_size_model_data %>%
  group_by(genotype, Isoline) %>%
  summarise(n = n())

overall <- sample_sizes %>%
  group_by(genotype) %>%
  summarise(n =sum(n)) %>%
  mutate(Isoline = "Across all isolines") %>%
  select(!! names(sample_sizes))

sample_sizes %>%
  bind_rows(overall) %>%
  kable(digits=3) %>% kable_styling()

Posterior model probabilities

The top model (with >99% probability) contains genotype and isoline (just the main effects). The model with an interaction, and the model lacking the isoline effect, had much lower probability. This indicated that SR/ST genotype affected the offspring sex ratio, and so did isoline, but there was little/no evidence for an interaction.

round(sort(SR_model_comparison, decreasing=T), 3) %>% 
  kable() %>% kable_styling()

Estimated offspring sex ratio by genotype

preds_all_isolines %>%
  map(~ posterior_summary(.x)) %>% do.call("rbind", .) %>% as_tibble() %>%
  mutate(Genotype = names(preds_all_isolines)) %>%
  select(Genotype, everything()) %>% 
  kable(digits = 3) %>% kable_styling()

Estimated offspring sex ratio by genotype and isoline

preds_by_isoline %>%
  map(~ posterior_summary(.x)) %>% do.call("rbind", .) %>% as_tibble() %>%
  mutate(Genotype = sex_ratio_data$genotype,
         Isoline = sex_ratio_data$isoline) %>%
  select(Genotype, Isoline, everything()) %>% 
  kable(digits = 3) %>% kable_styling()

Hypothesis testing

Here, we calculate the posterior difference in the mean sex ratio between each of the 3 genotypes:

get_p <- function(x) (100 - as.numeric(p_direction(x)))/100


data.frame(
  Comparison = c("SR/SR - ST/ST", "SR/ST - ST/ST", "SR/SR - SR/ST"),
  rbind(posterior_summary(preds_all_isolines$SRSR - preds_all_isolines$STST),
        posterior_summary(preds_all_isolines$SRST - preds_all_isolines$STST),
        posterior_summary(preds_all_isolines$SRSR - preds_all_isolines$SRST))
) %>% mutate(
  p = c(get_p(preds_all_isolines$SRSR - preds_all_isolines$STST),
        get_p(preds_all_isolines$SRST - preds_all_isolines$STST),
        get_p(preds_all_isolines$SRSR - preds_all_isolines$SRST))
) %>%
  kable(digits = 3) %>% kable_styling()

Figure showing sex ratios

SR_figure <- gather(preds_all_isolines) %>%
  ggplot(aes(key, value, fill = key)) + 
  geom_eye() + 
  geom_hline(yintercept = 50, linetype = 2) + 
  coord_cartesian(ylim = c(49, 65)) + 
  theme_bw() +
  theme(legend.position = "none") + 
  ylab("% daughters among adult offspring\n(posterior estimate)") + 
  xlab("Maternal genotype") 

SR_figure

Figure showing sex ratios - split by isoline

SR_figure2 <- gather(preds_by_isoline) %>%
  mutate(split = strsplit(key, split = "_"),
         Genotype = map_chr(split, ~.x[1]),
         isoline = paste(map_chr(split, ~.x[2]), "isoline")) %>%
  select(value, Genotype, isoline) %>%
  ggplot(aes(Genotype, value, fill = Genotype)) + 
  geom_eye() + 
  geom_hline(yintercept = 50, linetype = 2) + 
  facet_wrap(~isoline) + 
  coord_cartesian(ylim = c(48, 71.3)) + 
  theme_bw() +
  theme(legend.position = "none") + 
  ylab("% daughters among adult offspring\n(posterior estimate)") + 
  xlab("Maternal genotype")

SR_figure2

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    gridExtra_2.3      lmerTest_3.0-1    
 [4] lme4_1.1-17        Matrix_1.2-14      showtext_0.5-1    
 [7] showtextdb_2.0     sysfonts_0.7.2     RColorBrewer_1.1-2
[10] ggbeeswarm_0.6.0   kableExtra_0.9.0   bayestestR_0.2.2  
[13] brms_2.10.0        Rcpp_1.0.2         forcats_0.4.0     
[16] stringr_1.4.0      dplyr_0.8.3        purrr_0.3.2       
[19] readr_1.1.1        tidyr_0.8.2        tibble_2.1.3      
[22] ggplot2_3.1.0      tidyverse_1.2.1   

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