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)  
  )
  
  # 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"))
  
  model <- brm(model_formula,
               family = "hurdle_negbinomial",
               chains = 4, cores = 1, iter = 8000, 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")
  saveRDS(neat_R2(model), file = "output/R2_of_brms_model.rds")
} else model <- readRDS("output/brms_model.rds")

Graphically verify 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(model, type = "hist", nsamples = 11, binwidth = 5)

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

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 <- as.data.frame(summary(model)$random[[1]]) %>%
  rownames_to_column("Parameter") %>%
  mutate(p = NA,
         Parameter = c("sd(Block - Intercept)", "sd(Block - Hurdle intercept)"))

summary(model)$fixed %>% as.data.frame() %>% 
  rownames_to_column("Parameter") %>%
  left_join(bayesian_p_values, by = "Parameter") %>%
  rename(p = pd) %>% arrange(grepl("hu_", Parameter)) %>%
  rbind(random) %>%
  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, "*", ""),
         ` ` = replace(` `, is.na(` `), ""),
         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)
         ) %>% 
  kable() %>% kable_styling()

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.

# Define the data for prediction. We use the average values for continuous predictors like female age
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()

# Find the posterior medians and 95% CIs for each genotype/isoline combination
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)

# Get the posterior averaged across the isolines
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( # Summarise it, for Figure 1  
  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_fertile <- 
  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 * (1 - mean(.x))) %>% select(-genotype) %>% t() %>% as.data.frame() %>%
  rename(STST = V1, SRST = V2, SRSR = V3)
  
predicted_prop_fertile_all_iso <- data.frame(
  genotype = c("STST", "SRST", "SRSR"), Isoline = "All", 
  posterior_means_prop_fertile %>% lapply(posterior_summary) %>% do.call("rbind", .),
  facet = "C. % fertile females") 

# Concatenate the predictions into one data frame for plotting
preds_figure1 <- bind_rows(predicted_mean_all_iso,
                           predicted_mean_when_fertile_all_iso,
                           predicted_prop_fertile_all_iso) %>%
  mutate(genotype = factor(genotype, levels = c("STST", "SRST", "SRSR")))

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

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) 
  

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

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

figure_2