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)  
  )
  
  
  # 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)

  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(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")
} else{
  full_model <- readRDS("output/full_model.rds")
  genotype_only_model <- readRDS("output/genotype_only_model.rds")
  model_probabilities <- readRDS("output/model_comparison.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-5-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") + 
  geom_point(size = 3.1, pch = 21, colour = "black") + 
  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-9-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.