Statistical analysis

Selection treatment has affected activity and body size

There is a strong effect of selection treatment on activity in both sexes, and an effect of selection on female body size.

respiration %>%
  filter(CYCLE == "I") %>%
  mutate(`Body weight` = scale(BODY_WEIGHT),
         `Activity level` = scale(ACTIVITY)) %>%
  gather(trait, value, `Body weight`, `Activity level`) %>%
  ggplot(aes(SEX, value, colour = SELECTION)) + 
  stat_summary(position = position_dodge(0.4), size = 0.3) + 
  facet_wrap(~ trait) + 
  scale_color_brewer(palette = "Set1", direction = -1)

Scale the data

Here, I scale and centre the body mass (across all samples), and multiply VO2 and VCO2 by 1000 so that their units (and resulting regression coefficients) are close to those assumed by the brms default priors.

I did not scale and centre VO2 and VCO2, because I will soon relate them to each other via the respiratory quotient, RQ, so it makes sense to leave them in their original units. I also did not scale and centre activity level, because this variable is naturally bounded by zero and one, and so I can model it on its original scale using the beta distribution.

scaled_data <- respiration %>%
  mutate(VO2 = VO2 * 1000, 
         VCO2 = VCO2 * 1000,
         BODY_WEIGHT = as.numeric(scale(BODY_WEIGHT))) %>% 
  rename(BODYMASS = BODY_WEIGHT)

Activity level (one value per cycle, i.e. 3 measures on each ‘sample’ of individuals)

Formula: ACTIVITY ~ SELECTION * SEX + CYCLE + (1 | LINE) + (1 | SAMPLE)

This formula allows for effects on activity of sex and selection treatment (and their 2-way interaction), and for an effect of cycle (coded as a 3-level factor, allowing non-linear change across the 3 cycles). The random factors were added due to our repeated measures of replicate selection lines and samples (same for the following forrmulae).

Body mass

Formula: BODYMASS ~ SELECTION * SEX + (1 | LINE)

This formula allows for effects on activity of sex and selection treatment (and their 2-way interaction). Because there is only one measure of body mass for each sample of flies, we do not need to fit a sample-level random effect; also, this model is run on only a subset of the full dataset (one of the 3 cycles), since we would incur pseudo-replication if we used the full dataset.

VO2

Formula: VO2 ~ SELECTION * SEX * CYCLE + BODYMASS + ACTIVITY +
SELECTION:BODYMASS + SELECTION:ACTIVITY + SEX:BODYMASS + SEX:ACTIVITY + (1 | LINE) + (1 | SAMPLE)

This formula allows for effects on activity of sex, selection an cycle (and their 2- and 3-way interactions), and for sex- and selection treatment-specific effects of body mass and activity level.

VCO2 (as determined by the parameter RQ)

Formulae (2-part model, see vignette("brms_nonlinear")):

VCO2 ~ VO2 * (0.7 + 0.3 * inv_logit(RQ))

RQ ~ SELECTION * SEX * CYCLE + BODYMASS + ACTIVITY +
SELECTION:BODYMASS + SELECTION:ACTIVITY + SEX:BODYMASS + SEX:ACTIVITY + (1 | LINE) + (1 | SAMPLE)

VCO2 is assumed to depend on the value of VO2 from the same measurement, multiplied by RQ, a parameter that is constrained to vary between 0.7 and 1 (based on our prior knowledge of the chemistry of respiration) through the use of the inverse logit function. In turn, RQ is assumed to depend on the same set of predictors as for VO2.

Priors

To apply some mild regularisation and assist model convergence, we set a prior on all the fixed effect parameters of normal(0, 10).

Family

All response variables are assumed to follow a normal (Gaussian) distribution, except for activity level (beta); as we shall see, this turns out to be a reasonable approximation of the response variables’ true distributions.

Fit the brms model

inv_logit <- function(x) 1 / (1 + exp(-x))

# add a subsetting variable, so that we can estimate the effects of selection and sex on body size without having three redundant measures of body size (one per cycle). See ?brmsformula, section beginning "For multivariate models, subset may be used..."
scaled_data <- scaled_data %>%
      mutate(body_subset = CYCLE == "I")

if(!file.exists("output/brms_SEM.rds")){
  
  brms_formula <- 
    bf(VO2 ~ SELECTION * SEX * CYCLE + 
         BODYMASS + ACTIVITY +  
         SELECTION:BODYMASS + SELECTION:ACTIVITY + 
         SEX:BODYMASS + SEX:ACTIVITY + 
         (1 | LINE) + (1 | SAMPLE)) +
    bf(VCO2 ~ VO2 * (0.7 + 0.3 * inv_logit(RQ)), 
       RQ ~ SELECTION * SEX * CYCLE + 
         BODYMASS + ACTIVITY + 
         SELECTION:BODYMASS + SELECTION:ACTIVITY + 
         SEX:BODYMASS + SEX:ACTIVITY + 
         (1 | LINE) + (1 | SAMPLE),
       nl = TRUE) +
    bf(BODYMASS | subset(body_subset) ~ SELECTION * SEX + 
         (1 | LINE)) +
    bf(ACTIVITY ~ SELECTION * SEX + CYCLE + 
         (1 | LINE) + (1 | SAMPLE), family = "beta") +
    set_rescor(FALSE)
  
  brms_SEM <- brm(
    brms_formula,  
    data = scaled_data,
    iter = 10000, chains = 4, cores = 4,
    prior = prior(normal(0, 10), class = "b"),
    control = list(max_treedepth = 20, adapt_delta = 0.99)
  )
  
  saveRDS(brms_SEM, file = "output/brms_SEM.rds")
  bayesian_R2_values <- bayes_R2(brms_SEM)
  saveRDS(bayesian_R2_values,
          file = "output/bayesian_R2_values.rds")
  
} else {
  brms_SEM <- readRDS("output/brms_SEM.rds")
  bayesian_R2_values <-
    readRDS("output/bayesian_R2_values.rds")
}

Inspect the model output

Here is the complete output of summary() called on the model. Note that the model has converged ok (Rhat = 1), and that no parameters are under-sampled (shown by the ESS columns). Several parameters also differ significantly from zero (shown by their 95% credible intervals not overlapping zero). Note that the response variables are not all in the same units, so the magnitudes of their coefficients (“Estimate” column) are not directly comparable.

pvalues <- as.data.frame(p_direction(brms_SEM)) %>% 
  mutate(Parameter = str_remove_all(Parameter, "b_"),
         Parameter = str_replace_all(Parameter, "[.]", ":"),
         p = (100 - pd) / 100) %>%
  filter(!grepl("[.]1", Parameter))


fixef(brms_SEM) %>% 
  as.data.frame() %>%
  rownames_to_column("Parameter") %>%
  left_join(pvalues %>% select(-pd), by = "Parameter") %>%
  mutate(` ` = ifelse(p < 0.05, "\\*", ""),
         ` ` = replace(` `, p > 0.05 & p < 0.1, "~"),
         ` ` = replace(` `, p < 0.01, "**"), 
         ` ` = replace(` `, p < 0.001, "***")) %>%
  mutate(Response = map_chr(strsplit(Parameter, split = "_"), ~ .x[1]),
         Response = str_replace_all(Response, "BODYMASS", "Body mass"),
         Response = str_replace_all(Response, "ACTIVITY", "Activity"),
         Response = str_replace_all(Response, "VCO2", "RQ"),
         Parameter = str_replace_all(Parameter, "BODYMASS", "Body mass"),
         Parameter = str_replace_all(Parameter, "ACTIVITY", "Activity"),
         Parameter = str_remove_all(Parameter, ".+_"),
         Parameter = str_replace_all(Parameter, "SELECTIONPoly", "Polyandry"),
         Parameter = str_replace_all(Parameter, "CYCLEIII", "Cycle III"),
         Parameter = str_replace_all(Parameter, "CYCLEII", "Cycle II"),
         Parameter = str_replace_all(Parameter, "SEXM", "Male"),
         Parameter = str_replace_all(Parameter, ":", " x ")) %>%
  select(Response, Parameter, everything()) %>%
  mutate(Response = factor(Response, 
                           c("Activity", "Body mass", "VO2", "RQ"))) %>%
  arrange(Response) %>% select(-Response) %>%
  kable(digits = 3) %>% kable_styling() %>%
  group_rows("Activity level", 1, 6) %>%
  group_rows("Body mass", 7, 10) %>%
  group_rows("VO2", 11, 28) %>%
  group_rows("Respiratory quotient (RQ)", 29, 46)

Inspect the model R^2 values

These are all pretty high, which doesn’t tell us much, but it’s included here for completeness.

bayesian_R2_values %>%
  as.data.frame() %>%
  rownames_to_column("Response variable") %>%
  mutate(`Response variable` = gsub("R2", "", `Response variable`)) %>%
  kable(digits=2) %>% kable_styling()

Plot posterior predictive checks

Here, we draw 10 samples from the posterior of the model, and apply them to the original predictor variables, to check that the predicted values resemble the real data (they should, if the model is a sensible approximation of the real processes that generated the data). The plot looks good, mostly, though the fit is a little off for ACTIVITY because this variable is bounded at zero (prior to being scaled). Rather than adding a more complex model, we elected to settle for the Gaussian approximation.

my_check <- function(r) {
  pp_check(brms_SEM, resp = r) + 
    theme(legend.position = "none") + 
    labs(subtitle = r)
}

grid.arrange(
  my_check("VO2"),
  my_check("VCO2"),
  my_check("ACTIVITY"),
  my_check("BODYMASS")
)

Extract the posterior estimates of the means

These are used for plotting the range of means that is supported by the data, given our priors. The posterior estimates show the mean of each group, accounting for all the random effects (i.e. the design of the experiment), the covariance structure of the response variables, etc.

We will also use these posteriors for hypothesis testing, e.g. to see if the mean body size of the polyandry flies differs from that of monogamy flies, by subtracting one posterior from the other to get the posterior estimate of the difference in means. If most (e.g. >95%) of this posterior difference lies on one side of zero, the two means may be considered ‘significantly different’ as conventionally defined. The magnitude of the difference in means is also an intuitive measure of effect size, and the posterior gives a sense of how precisely we have estimated effect size.

new <- scaled_data %>%
  select(SELECTION, SEX, CYCLE) %>%
  distinct() %>%
  mutate(body_subset = TRUE)

# Get the posterior estimate of body weight and activity, 
# for the 4 combinations of sex and monogamy/polandry treatment
bodymass <- brms_SEM %>% 
  fitted(new, re_formula = NA, 
         resp = "BODYMASS", summary = FALSE)
activity <- brms_SEM %>% 
  fitted(new, re_formula = NA, 
         resp = "ACTIVITY", summary = FALSE)

if(!file.exists("data/posterior_predictions.rds")){
  
  # Calculate the posterior VO2 and RQ for the levels of 'new', either by assuming bodysize and acitivity are at their global means, or by using the means we have estimated for each combination of selection, sex, and cycle
  get_post_response <- function(body, act){
    
    n_combos <- nrow(new)
    if(!body) bodyweight <- matrix(0, ncol = n_combos, nrow = 1)
    if(!act)  activity <- matrix(0, ncol = n_combos, nrow = 1)
    
    get_diagonal <- function(mat){
      if(ncol(mat) == nrow(mat)) return(diag(mat))
      mat
    }
    
    VO2_posterior <- map(
      1:n_combos, 
      ~ brms_SEM %>% 
        fitted(data.frame(
          SELECTION = new$SELECTION[.x],
          SEX = new$SEX[.x],
          CYCLE = new$CYCLE[.x],
          BODYMASS = bodymass[, .x], 
          ACTIVITY = activity[, .x]), 
          re_formula = NA, resp = "VO2", 
          summary = FALSE) %>%
        get_diagonal()) %>%
      do.call("cbind", .)
    
    RQ_posterior <- map(
      1:n_combos, 
      ~ brms_SEM %>% 
        fitted(data.frame(
          SELECTION = new$SELECTION[.x],
          SEX = new$SEX[.x],
          BODYMASS = bodymass[, .x], 
          ACTIVITY = activity[, .x],
          CYCLE = new$CYCLE[.x],
          VO2 = VO2_posterior[, .x]), 
          re_formula = NA, resp = "VCO2", 
          nlpar = "RQ", summary = FALSE) %>%
        get_diagonal()) %>%
      do.call("cbind", .)
    
    n_draws <- nrow(VO2_posterior)
    
    posterior <- tibble(
      value = c(VO2_posterior),
      resp = "VO2",
      selection = rep(new$SELECTION, each = n_draws),
      sex = rep(new$SEX, each = n_draws),
      cycle = rep(new$CYCLE, each = n_draws)
    ) %>% bind_rows(
      tibble(
        value = c(RQ_posterior),
        resp = "RQ",
        selection = rep(new$SELECTION, each = n_draws),
        sex = rep(new$SEX, each = n_draws),
        cycle = rep(new$CYCLE, each = n_draws)
      )
    ) %>% mutate(
      value = ifelse(resp == "RQ", 
                     0.7 + 0.3 * inv_logit(value), 
                     value))
  }
  
  posterior_predictions <- list(
    direct_path = get_post_response(body = F, act = F),
    direct_path_plus_bodyweight = get_post_response(body = T, act = F),
    direct_path_plus_activity = get_post_response(body = F, act = T),
    all_paths = get_post_response(body = T, act = T))
  
  posterior_predictions <- lapply(
    names(posterior_predictions), function(i) {
      posterior_predictions[[i]] %>% mutate(path = i)
    }) %>% bind_rows() %>%
    mutate(path = factor(path, unique(path)))
  
  posterior_predictions <- posterior_predictions %>%
    mutate(sex = ifelse(sex == "M", "Male", "Female"),
           sex = factor(sex, c("Male", "Female")),
           path = replace(as.character(path), 
                          path == "direct_path", 
                          "1. Direct effect"),
           path = replace(path, 
                          path == "direct_path_plus_bodyweight", 
                          "2. Direct + body mass"),
           path = replace(path, 
                          path == "direct_path_plus_activity", 
                          "3. Direct + activity level"),
           path = replace(path, 
                          path == "all_paths", 
                          "4. Direct + both mediators")) 
  
  posterior_predictions %>% 
    saveRDS("data/posterior_predictions.rds")
} else {
  posterior_predictions <- 
    readRDS("data/posterior_predictions.rds")
}

Plot the posterior means for VO2

The plot show the model’s best estimates of the mean, with the shading showing the posterior distribution of these guesses. The four rows show the effect of: - selection on VO2, independently of both mediator variables - selection on VO2, independently of activity level but incorporating the mediating effect of body mass - selection on VO2, independently of body weight but incorporating the mediating effect of activity level - selection on VO2, including the effects of both mediator variables

posterior_predictions %>% 
  filter(resp == "VO2") %>%
  ggplot(aes(x = cycle, y = scale(value), fill = selection)) +
  scale_fill_brewer(palette = "Set1", direction = -1, name = NULL) +
  geom_eye (alpha=0.6, .width = c(0.88), position = pd) +
  coord_cartesian(ylim = c(-3.5, 2.8)) +
  facet_grid(path~sex) + 
  theme_bw() + 
  theme(legend.position = "top",
        panel.grid.major.x = element_blank()) +   
  xlab("Cycle") + 
  ylab("Oxygen consumption (Posterior estimate of mean VO2)") 

Plot the posterior means for RQ

posterior_predictions %>%
  filter(resp == "RQ") %>%
  ggplot(aes(x = cycle, y = value, fill = selection)) +
  scale_fill_brewer(palette = "Set1", direction = -1) +
  geom_eye (alpha=0.6, .width = c(0.88), position = pd) +
  facet_grid(path~sex) + 
  theme_bw() + 
  theme(legend.position = "top",
        panel.grid.major.x = element_blank()) +   
  xlab("Cycle") + 
  ylab("RQ (Posterior estimate of group means)")