Metabolite composition
Last updated: 2021-10-06
Checks: 7 0
Knit directory: exp_evol_respiration/
This reproducible R Markdown analysis was created with workflowr (version 1.6.2). The Checks tab describes the reproducibility checks that were applied when the results were created. The Past versions tab lists the development history.
Great! Since the R Markdown file has been committed to the Git repository, you know the exact version of the code that produced these results.
Great job! The global environment was empty. Objects defined in the global environment can affect the analysis in your R Markdown file in unknown ways. For reproduciblity it’s best to always run the code in an empty environment.
The command set.seed(20190703) was run prior to running the code in the R Markdown file. Setting a seed ensures that any results that rely on randomness, e.g. subsampling or permutations, are reproducible.
Great job! Recording the operating system, R version, and package versions is critical for reproducibility.
Nice! There were no cached chunks for this analysis, so you can be confident that you successfully produced the results during this run.
Great job! Using relative paths to the files within your workflowr project makes it easier to run your code on other machines.
Great! You are using Git for version control. Tracking code development and connecting the code version to the results is critical for reproducibility.
The results in this page were generated with repository version 7095d8f. See the Past versions tab to see a history of the changes made to the R Markdown and HTML files.
Note that you need to be careful to ensure that all relevant files for the analysis have been committed to Git prior to generating the results (you can use wflow_publish or wflow_git_commit). workflowr only checks the R Markdown file, but you know if there are other scripts or data files that it depends on. Below is the status of the Git repository when the results were generated:
Ignored files:
Ignored: .DS_Store
Ignored: .Rhistory
Ignored: .Rproj.user/
Ignored: data/.DS_Store
Ignored: figures/.DS_Store
Ignored: output/.DS_Store
Untracked files:
Untracked: data/1.eclosion_wide.csv
Untracked: data/2.DesRes.csv
Untracked: data/2.StarvRes.csv
Untracked: data/3.metabolic_rates.csv
Untracked: data/4.metabolite_data.csv
Untracked: figures/desiccation.pdf
Untracked: figures/eclosion.pdf
Untracked: figures/eclosion_interaction_plot.pdf
Untracked: figures/metaboliteDAG.pdf
Untracked: figures/metabolite_plotCONTROLLED.pdf
Untracked: figures/metabolite_plotCONTROLLED_nolines.pdf
Untracked: figures/metabolite_plot_nolines.pdf
Untracked: figures/metabolites_supp.pdf
Untracked: figures/respirationDAG.pdf
Untracked: figures/respiration_figure_cycleI.pdf
Untracked: figures/respiration_figure_cycleI_inkscape.pdf
Untracked: figures/respiration_figure_noRQ.pdf
Untracked: figures/starvation.pdf
Untracked: output/brms_SEM_respiration_noslope.rds
Untracked: output/brms_metabolite_SEM_noslope.rds
Untracked: output/cox_brms_noslope.rds
Untracked: output/des_brm_noslope.rds
Untracked: output/metabolite_SEM_table.csv
Untracked: output/respirometry_SEM_table.csv
Untracked: output/sta_brm_noslope.rds
Untracked: output/stress_medians.csv
Untracked: output/wing_brms_noslope.rds
Unstaged changes:
Modified: analysis/respirometry.Rmd
Deleted: data/1.eclosion_times.csv
Deleted: data/2.metabolic_rates.csv
Deleted: data/3.DesRes.csv
Deleted: data/3.StarvRes.csv
Deleted: data/3.metabolite_data.csv
Modified: figures/metabolite_interaction_plot.pdf
Modified: figures/metabolite_pairs_plot.pdf
Modified: figures/metabolite_plot.pdf
Modified: figures/respiration_figure.pdf
Modified: figures/respiration_pairs_plot.pdf
Modified: output/cox_brms.rds
Deleted: output/coxmod.rds
Deleted: output/coxmod_dropINT.rds
Deleted: output/coxmod_dropSEX.rds
Deleted: output/coxmod_dropTRT.rds
Note that any generated files, e.g. HTML, png, CSS, etc., are not included in this status report because it is ok for generated content to have uncommitted changes.
These are the previous versions of the repository in which changes were made to the R Markdown (analysis/metabolites.Rmd) and HTML (docs/metabolites.html) files. If you’ve configured a remote Git repository (see ?wflow_git_remote), click on the hyperlinks in the table below to view the files as they were in that past version.
| File | Version | Author | Date | Message |
|---|---|---|---|---|
| Rmd | 7095d8f | MartinGarlovsky | 2021-10-06 | wflow_publish(“analysis/metabolites.Rmd”) |
| html | aa3e11e | MartinGarlovsky | 2021-10-05 | Build site. |
| Rmd | a823c3e | MartinGarlovsky | 2021-10-05 | wflow_publish(c(“analysis/metabolites.Rmd”)) |
| html | 21567bb | MartinGarlovsky | 2021-03-12 | Build site. |
| Rmd | f7b49ff | MartinGarlovsky | 2021-03-12 | wflow_publish(c(“analysis/juvenile_development.Rmd”, “analysis/resistance.Rmd”, |
| html | ffb09dd | MartinGarlovsky | 2021-02-08 | Build site. |
| Rmd | 1b92046 | MartinGarlovsky | 2021-02-08 | wflow_publish(c(“analysis/juvenile_development.Rmd”, “analysis/resistance.Rmd”, |
| html | 41d232f | lukeholman | 2020-12-18 | Build site. |
| Rmd | 99fd9d4 | lukeholman | 2020-12-18 | new menu |
| html | 0d5bcc9 | lukeholman | 2020-12-18 | Build site. |
| html | 989e86f | lukeholman | 2020-12-18 | Build site. |
| html | a82b6fb | lukeholman | 2020-12-18 | Build site. |
| Rmd | c058bc3 | lukeholman | 2020-12-18 | new menu |
| html | b1f42b1 | lukeholman | 2020-12-18 | Build site. |
| Rmd | bcaef70 | lukeholman | 2020-12-18 | Neater figures |
| html | 034431b | lukeholman | 2020-12-18 | Build site. |
| Rmd | a3b516e | lukeholman | 2020-12-18 | Neater figures |
| html | b6cf554 | lukeholman | 2020-12-11 | Build site. |
| Rmd | dcb0a67 | lukeholman | 2020-12-11 | Tweaks |
| html | a88f037 | lukeholman | 2020-12-11 | Build site. |
| Rmd | 9588f44 | lukeholman | 2020-12-11 | Tweaks |
| html | bb96acf | lukeholman | 2020-12-11 | Build site. |
| html | e5c580f | lukeholman | 2020-12-11 | Build site. |
| Rmd | 328e117 | lukeholman | 2020-12-11 | Tweaks |
| html | 49e7792 | lukeholman | 2020-12-11 | Build site. |
| Rmd | 9d65fc7 | lukeholman | 2020-12-11 | Tweaks |
| html | 2911fb9 | lukeholman | 2020-12-10 | Build site. |
| Rmd | 184b0a4 | lukeholman | 2020-12-10 | Tweaks |
| html | a361e42 | lukeholman | 2020-12-10 | Build site. |
| Rmd | 98c43b5 | lukeholman | 2020-12-10 | Tweaks |
| html | 7fca240 | lukeholman | 2020-12-10 | Build site. |
| Rmd | 125d4a2 | lukeholman | 2020-12-10 | Tweaks |
| html | 31ed22a | lukeholman | 2020-12-10 | Build site. |
| Rmd | 8464a7f | lukeholman | 2020-12-10 | Tweaks |
| html | 901053c | lukeholman | 2020-12-10 | Build site. |
| Rmd | 904af31 | lukeholman | 2020-12-10 | Tweaks |
| html | f7c88a2 | lukeholman | 2020-12-10 | Build site. |
| Rmd | 68780f6 | lukeholman | 2020-12-10 | Tweaks |
| html | deb7183 | lukeholman | 2020-12-09 | Build site. |
| Rmd | 720eb6d | lukeholman | 2020-12-09 | Tweaks |
| html | b731971 | lukeholman | 2020-12-09 | Build site. |
| Rmd | 398d963 | lukeholman | 2020-12-09 | Tweaks |
| html | b449eb3 | lukeholman | 2020-12-09 | Build site. |
| Rmd | 15f3c92 | lukeholman | 2020-12-09 | Tweaks |
| html | 43cc270 | lukeholman | 2020-12-09 | Build site. |
| Rmd | 2642c27 | lukeholman | 2020-12-09 | More work |
| html | 4f5ee28 | lukeholman | 2020-12-04 | Build site. |
| Rmd | d441b69 | lukeholman | 2020-12-04 | Luke metabolites analysis |
| Rmd | c8feb2d | lukeholman | 2020-11-30 | Same page with Martin |
| html | 3fdbcb2 | lukeholman | 2020-11-30 | Tweaks Nov 2020 |
Load packages
# it was slightly harder to install the showtext package. On Mac, I did this:
# installed 'homebrew' using Terminal: ruby -e "$(curl -fsSL https://raw.githubusercontent.com/Homebrew/install/master/install)"
# installed 'libpng' using Terminal: brew install libpng
# installed 'showtext' in R using: devtools::install_github("yixuan/showtext")
library(tidyverse)
library(GGally)
library(gridExtra)
library(ggridges)
library(brms)
library(tidybayes)
library(DT)
library(kableExtra)
library(knitrhooks) # install with devtools::install_github("nathaneastwood/knitrhooks")
library(showtext)
output_max_height() # a knitrhook option
# set up nice font for figure
nice_font <- "Lora"
font_add_google(name = nice_font, family = nice_font, regular.wt = 400, bold.wt = 700)
showtext_auto()
options(stringsAsFactors = FALSE)Load metabolite composition data
This analysis set out to test whether sexual selection treatment had an effect on macro-metabolite composition of flies. We measured fresh and dry fly weight in milligrams, plus the concentrations of five metabolites. These are:
Lipid_conc(i.e. the weight of the hexane fraction, divided by the full dry weight),Carbohydrate_conc(i.e. the weight of the aqueous fraction, divided by the full dry weight),Protein_conc(i.e. \(\mu\)g of protein per milligram as measured by the bicinchoninic acid protein assay),Glycogen_conc(i.e. \(\mu\)g of glycogen per milligram as measured by the hexokinase assay), andChitin_conc(estimated as the difference between the initial and final dry weights)
We expect body weight to vary between the sexes and potentially between treatments. In turn, we expect body weight might co-vary with the five metabolite concentrations, e.g. because flies might become heavier by sequestering proportionally more of particular metabolites. There might also be weight-independent effects of sex and selection treatment on metabolite composition.
metabolites <- read_csv('data/4.metabolite_data.csv') %>%
mutate(sex = ifelse(sex == "m", "Male", "Female"),
line = paste(treatment, line, sep = ""),
treatment = ifelse(treatment == "M", "Monogamy", "Polyandry")) %>%
# log transform glycogen since it shows a long tail (others are reasonably normal-looking)
mutate(Glycogen_ug_mg = log(Glycogen_ug_mg)) %>%
# There was a technical error with flies collected on day 1,
# so they are excluded from the whole paper. All the measurements analysed are of 3d-old flies
filter(time == '2') %>%
select(-time)
scaled_metabolites <- metabolites %>%
# Find proportional metabolites as a proportion of total dry weight
mutate(
Dry_weight = dwt_mg,
Lipid_conc = Hex_frac / Dry_weight,
Carbohydrate_conc = Aq_frac / Dry_weight,
Protein_conc = Protein_ug_mg,
Glycogen_conc = Glycogen_ug_mg,
Chitin_conc = Chitin_mg_mg) %>%
select(sex, treatment, line, Dry_weight, ends_with("conc")) %>%
mutate_at(vars(ends_with("conc")), ~ as.numeric(scale(.x))) %>%
mutate(Dry_weight = as.numeric(scale(Dry_weight))) %>%
mutate(sextreat = paste(sex, treatment),
sextreat = replace(sextreat, sextreat == "Male Monogamy", "M males"),
sextreat = replace(sextreat, sextreat == "Male Polyandry", "E males"),
sextreat = replace(sextreat, sextreat == "Female Monogamy", "M females"),
sextreat = replace(sextreat, sextreat == "Female Polyandry", "E females"),
sextreat = factor(sextreat, c("M males", "E males", "M females", "E females")))Inspect the raw data
Raw numbers
All variables are shown in standard units (i.e. mean = 0, SD = 1).
my_data_table <- function(df){
datatable(
df, rownames=FALSE,
autoHideNavigation = TRUE,
extensions = c("Scroller", "Buttons"),
options = list(
dom = 'Bfrtip',
deferRender=TRUE,
scrollX=TRUE, scrollY=400,
scrollCollapse=TRUE,
buttons =
list('csv', list(
extend = 'pdf',
pageSize = 'A4',
orientation = 'landscape',
filename = 'Dpseudo_metabolites')),
pageLength = 50
)
)
}
scaled_metabolites %>%
select(-sextreat) %>%
mutate_if(is.numeric, ~ format(round(.x, 3), nsmall = 3)) %>%
my_data_table()Simple plots
The following plot shows how each metabolite varies between sexes and treatments, and how the concentration of each metabolite co-varies with dry weight across individuals.
levels <- c("Carbohydrate", "Chitin", "Glycogen", "Lipid", "Protein", "Dry weight")
cols <- c("M females" = "pink",
"E females" = "red",
"M males" = "skyblue",
"E males" = "blue")
grid.arrange(
scaled_metabolites %>%
rename_all(~ str_remove_all(.x, "_conc")) %>%
rename(`Dry weight` = Dry_weight) %>%
mutate(sex = factor(sex, c("Male", "Female"))) %>%
reshape2::melt(id.vars = c('sex', 'treatment', 'sextreat', 'line')) %>%
mutate(variable = factor(variable, levels)) %>%
ggplot(aes(x = sex, y = value, fill = sextreat)) +
geom_hline(yintercept = 0, linetype = 2) +
geom_boxplot() +
facet_grid( ~ variable) +
theme_bw() +
xlab("Sex") + ylab("Concentration") +
theme(legend.position = 'top',
text = element_text(family = nice_font)) +
scale_fill_manual(values = cols, name = ""),
arrangeGrob(
scaled_metabolites %>%
rename_all(~ str_remove_all(.x, "_conc")) %>%
reshape2::melt(id.vars = c('sex', 'treatment', 'sextreat', 'line', 'Dry_weight')) %>%
mutate(variable = factor(variable, levels)) %>%
ggplot(aes(x = Dry_weight, y = value, colour = sextreat, fill = sextreat)) +
geom_smooth(method = 'lm', se = TRUE, aes(colour = NULL, fill = NULL), colour = "grey20", size = .4) +
geom_point(pch = 21, colour = "grey20") +
facet_grid( ~ variable) +
theme_bw() +
xlab("Dry weight") + ylab("Concentration") +
theme(legend.position = 'none',
text = element_text(family = nice_font),) +
scale_colour_manual(values = cols, name = "") +
scale_fill_manual(values = cols, name = ""),
scaled_metabolites %>%
rename_all(~ str_remove_all(.x, "_conc")) %>%
reshape2::melt(id.vars = c('sex', 'treatment', 'sextreat', 'line', 'Dry_weight')) %>%
mutate(variable = factor(variable, levels)) %>%
ggplot(aes(x = Dry_weight, y = value, colour = sextreat, fill = sextreat)) +
theme_void() + ylab(NULL), nrow = 1, widths = c(0.84, 0.16)),
heights = c(0.55, 0.45)
)
Plot of correlations between variables
Some of the metabolites, especially lipid concentration, are correlated with dry weight. There is also a large difference in dry weight between sexes (and treatments, to a lesser extent), and sex and treatment effects are evident for some of the metabolites in the raw data. Some of the metabolites are weakly correlated with other metabolites, e.g. lipid and glycogen concentration.
modified_densityDiag <- function(data, mapping, ...) {
ggally_densityDiag(data, mapping, colour = "grey10", ...) +
scale_fill_manual(values = cols) +
scale_x_continuous(guide = guide_axis(check.overlap = TRUE))
}
modified_points <- function(data, mapping, ...) {
ggally_points(data, mapping, pch = 21, colour = "grey10", ...) +
scale_fill_manual(values = cols) +
scale_x_continuous(guide = guide_axis(check.overlap = TRUE))
}
modified_facetdensity <- function(data, mapping, ...) {
ggally_facetdensity(data, mapping, ...) +
scale_colour_manual(values = cols)
}
modified_box_no_facet <- function(data, mapping, ...) {
ggally_box_no_facet(data, mapping, colour = "grey10", ...) +
scale_fill_manual(values = cols)
}
metabolite_pairs_plot <- scaled_metabolites %>%
arrange(sex, treatment) %>%
select(-line, -sex, -treatment) %>%
rename(`Sex and treatment` = sextreat) %>%
rename_all(~ str_replace_all(.x, "_", " ")) %>%
ggpairs(aes(colour = `Sex and treatment`, fill = `Sex and treatment`),
diag = list(continuous = wrap(modified_densityDiag, alpha = 0.7),
discrete = wrap("blank")),
lower = list(continuous = wrap(modified_points, alpha = 0.7, size = 1.1),
discrete = wrap("blank"),
combo = wrap(modified_box_no_facet, alpha = 0.7)),
upper = list(continuous = wrap(modified_points, alpha = 0.7, size = 1.1),
discrete = wrap("blank"),
combo = wrap(modified_box_no_facet, alpha = 0.7, size = 0.5)))
#metabolite_pairs_plot %>% ggsave(filename = "figures/metabolite_pairs_plot.pdf", height = 10, width = 10)
metabolite_pairs_plot
Mean dry weight
se <- function(x) sd(x) / sqrt(length(x))
metabolites %>%
group_by(sex, treatment) %>%
summarise(mean_dwt = mean(dwt_mg),
SE = se(dwt_mg),
n = n()) %>%
kable(digits = 3) %>% kable_styling(full_width = FALSE)| sex | treatment | mean_dwt | SE | n |
|---|---|---|---|---|
| Female | Monogamy | 0.562 | 0.019 | 12 |
| Female | Polyandry | 0.644 | 0.017 | 12 |
| Male | Monogamy | 0.330 | 0.009 | 12 |
| Male | Polyandry | 0.353 | 0.009 | 12 |
Directed acyclic graph (DAG)
This directed acyclic graph (DAG) illustrates the causal pathways that we observed between the experimental or measured variables (square boxes) and latent variables (ovals). We hypothesise that sex and mating system potentially influence dry weight as well as the metabolite composition (which we assessed by estimating the concentrations of carbohydrates, chitin, glycogen, lipids and protein). Additionally, dry weight is likely correlated with metabolite composition, and so dry weight acts as a ‘mediator variable’ between metabolite composition, and sex and treatment. The structural equation model below is built with this DAG in mind.
DiagrammeR::grViz('digraph {
graph [layout = dot, rankdir = LR]
# define the global styles of the nodes. We can override these in box if we wish
node [shape = rectangle, style = filled, fillcolor = Linen]
"Metabolite\ncomposition" [shape = oval, fillcolor = Beige]
# edge definitions with the node IDs
"Mating system\ntreatment (E vs M)" -> {"Dry weight"}
"Mating system\ntreatment (E vs M)" -> {"Metabolite\ncomposition"}
"Sex\n(Female vs Male)" -> {"Dry weight"} -> {"Metabolite\ncomposition"}
"Sex\n(Female vs Male)" -> {"Metabolite\ncomposition"}
{"Metabolite\ncomposition"} -> "Carbohydrates"
{"Metabolite\ncomposition"} -> "Chitin"
{"Metabolite\ncomposition"} -> "Glycogen"
{"Metabolite\ncomposition"} -> "Lipids"
{"Metabolite\ncomposition"} -> "Protein"
}')Fit brms structural equation model
Here we fit a model of the five metabolites, which includes dry body weight as a mediator variable. That is, our model estimates the effect of treatment, sex and line (and all the 2- and 3-way interactions) on dry weight, and then estimates the effect of those some predictors (plus dry weight) on the five metabolites. The model assumes that although the different sexes, treatment groups, and lines may differ in their dry weight, the relationship between dry weight and the metabolites does not vary by sex/treatment/line. This assumption was made to constrain the number of parameters in the model, and to reflect out prior beliefs about allometric scaling of metabolites.
Define Priors
We set fairly tight Normal priors on all fixed effect parameters, which ‘regularises’ the estimates towards zero – this is conservative (because it ensures that a stronger signal in the data is needed to produce a given posterior effect size estimate), and it also helps the model to converge. Similarly, we set a somewhat conservative half-cauchy prior (mean 0, scale 0.01) on the random effects for line (i.e. we consider large differences between lines – in terms of means and treatment effects – to be possible but improbable). We leave all other priors at the defaults used by brms. Note that the Normal priors are slightly wider in the model of dry weight, because we expect larger effect sizes of sex and treatment on dry weight than on the metabolite composition.
prior1 <- c(set_prior("normal(0, 0.5)", class = "b", resp = 'Lipid'),
set_prior("normal(0, 0.5)", class = "b", resp = 'Carbohydrate'),
set_prior("normal(0, 0.5)", class = "b", resp = 'Protein'),
set_prior("normal(0, 0.5)", class = "b", resp = 'Glycogen'),
set_prior("normal(0, 0.5)", class = "b", resp = 'Chitin'),
set_prior("normal(0, 1)", class = "b", resp = 'Dryweight'),
set_prior("cauchy(0, 0.01)", class = "sd", resp = 'Lipid', group = "line"),
set_prior("cauchy(0, 0.01)", class = "sd", resp = 'Carbohydrate', group = "line"),
set_prior("cauchy(0, 0.01)", class = "sd", resp = 'Protein', group = "line"),
set_prior("cauchy(0, 0.01)", class = "sd", resp = 'Glycogen', group = "line"),
set_prior("cauchy(0, 0.01)", class = "sd", resp = 'Chitin', group = "line"),
set_prior("cauchy(0, 0.01)", class = "sd", resp = 'Dryweight', group = "line"))
prior1 prior class coef group resp dpar nlpar bound source
normal(0, 0.5) b Lipid user
normal(0, 0.5) b Carbohydrate user
normal(0, 0.5) b Protein user
normal(0, 0.5) b Glycogen user
normal(0, 0.5) b Chitin user
normal(0, 1) b Dryweight user
cauchy(0, 0.01) sd line Lipid user
cauchy(0, 0.01) sd line Carbohydrate user
cauchy(0, 0.01) sd line Protein user
cauchy(0, 0.01) sd line Glycogen user
cauchy(0, 0.01) sd line Chitin user
cauchy(0, 0.01) sd line Dryweight user
Define the SEM’s sub-models
The fixed effects formula is sex * treatment + Dryweight (or sex * treatment in the case of the model of dry weight). The random effects part of the formula indicates that the 8 independent selection lines may differ in their means, and that the treatment effect may vary in sign/magnitude between lines. The notation | p | means that the model estimates the correlations in line effects (both slopes and intercepts) between the 6 response variables. Finally, the notation set_rescor(TRUE) means that the model should estimate the residual correlations between the response variables.
brms_formula <-
# Sub-models of the 5 metabolites
bf(mvbind(Lipid, Carbohydrate, Protein, Glycogen, Chitin) ~
sex*treatment + Dryweight + (treatment | p | line)) +
# dry weight sub-model
bf(Dryweight ~ sex*treatment + (treatment | p | line)) +
# Allow for (and estimate) covariance between the residuals of the difference response variables
set_rescor(TRUE)
brms_formulaLipid ~ sex * treatment + Dryweight + (treatment | p | line) Carbohydrate ~ sex * treatment + Dryweight + (treatment | p | line) Protein ~ sex * treatment + Dryweight + (treatment | p | line) Glycogen ~ sex * treatment + Dryweight + (treatment | p | line) Chitin ~ sex * treatment + Dryweight + (treatment | p | line) Dryweight ~ sex * treatment + (treatment | p | line)
Running the model
The model is run over 4 chains with 5000 iterations each (with the first 2500 discarded as burn-in), for a total of 2500*4 = 10,000 posterior samples.
if(!file.exists("output/brms_metabolite_SEM.rds")){
brms_metabolite_SEM <- brm(
brms_formula,
data = scaled_metabolites %>% # brms does not like underscores in variable names
rename(Dryweight = Dry_weight) %>%
rename_all(~ gsub("_conc", "", .x)),
iter = 5000, chains = 4, cores = 1,
prior = prior1,
control = list(max_treedepth = 20,
adapt_delta = 0.99)
)
#saveRDS(brms_metabolite_SEM, "output/brms_metabolite_SEM_noslope.rds") # save with no random slope term
saveRDS(brms_metabolite_SEM, "output/brms_metabolite_SEM.rds")
} else {
brms_metabolite_SEM <- readRDS('output/brms_metabolite_SEM_noslope.rds')
}Posterior predictive check of model fit
The plot below shows that the fitted model is able to produce posterior predictions that have a similar distribution to the original data, for each of the response variables, which is a necessary condition for the model to be used for statistical inference.
grid.arrange(
pp_check(brms_metabolite_SEM, resp = "Dryweight") +
ggtitle("Dry weight") + theme(legend.position = "none"),
pp_check(brms_metabolite_SEM, resp = "Lipid") +
ggtitle("Lipid") + theme(legend.position = "none"),
pp_check(brms_metabolite_SEM, resp = "Carbohydrate") +
ggtitle("Carbohydrate") + theme(legend.position = "none"),
pp_check(brms_metabolite_SEM, resp = "Protein") +
ggtitle("Protein") + theme(legend.position = "none"),
pp_check(brms_metabolite_SEM, resp = "Glycogen") +
ggtitle("Glycogen") + theme(legend.position = "none"),
pp_check(brms_metabolite_SEM, resp = "Chitin") +
ggtitle("Chitin") + theme(legend.position = "none"),
nrow = 2
)
Table of model parameter estimates
Formatted table
This tables shows the fixed effects estimates for treatment, sex, their interaction, as well as the slope associated with dry weight (where relevant), for each of the six response variables. The p column shows 1 - minus the “probability of direction”, i.e. the posterior probability that the reported sign of the estimate is correct given the data and the prior; subtracting this value from one gives a Bayesian equivalent of a one-sided p-value. For brevity, we have omitted all the parameter estimates involving the predictor variable line, as well as the estimates of residual (co)variance. Click the next tab to see a complete summary of the model and its output.
vars <- c("Lipid", "Carbohydrate", "Glycogen", "Protein", "Chitin")
tests <- c('_Dryweight', '_sexMale',
'_sexMale:treatmentPolyandry',
'_treatmentPolyandry')
hypSEM <- data.frame(expand_grid(vars, tests) %>%
mutate(est = NA,
err = NA,
lwr = NA,
upr = NA) %>%
# bind body weight on the end
rbind(data.frame(
vars = rep('Dryweight', 3),
tests = c('_sexMale',
'_treatmentPolyandry',
'_sexMale:treatmentPolyandry'),
est = NA,
err = NA,
lwr = NA,
upr = NA)))
for(i in 1:nrow(hypSEM)) {
result = hypothesis(brms_metabolite_SEM,
paste0(hypSEM[i, 1], hypSEM[i, 2], ' = 0'))$hypothesis
hypSEM[i, 3] = round(result$Estimate, 3)
hypSEM[i, 4] = round(result$Est.Error, 3)
hypSEM[i, 5] = round(result$CI.Lower, 3)
hypSEM[i, 6] = round(result$CI.Upper, 3)
}
pvals <- bayestestR::p_direction(brms_metabolite_SEM) %>%
as.data.frame() %>%
mutate(vars = map_chr(str_split(Parameter, "_"), ~ .x[2]),
tests = map_chr(str_split(Parameter, "_"), ~ .x[3]),
tests = str_c("_", str_remove_all(tests, "[.]")),
tests = replace(tests, tests == "_sexMaletreatmentPolyandry", "_sexMale:treatmentPolyandry")) %>%
filter(!str_detect(tests, "line")) %>%
mutate(p_val = 1 - pd, star = ifelse(p_val < 0.05, "\\*", "")) %>%
select(vars, tests, p_val, star)
hypSEM <- hypSEM %>% left_join(pvals, by = c("vars", "tests"))
hypSEM %>%
mutate(Parameter = c(rep(c('Dry weight', 'Sex (M)',
'Sex (M) x Treatment (E)',
'Treatment (E)'), 5),
'Sex (M)', 'Treatment (E)', 'Sex (M) x Treatment (E)')) %>%
mutate(Parameter = factor(Parameter, c("Dry weight", "Sex (M)", "Treatment (E)", "Sex (M) x Treatment (E)")),
vars = factor(vars, c("Carbohydrate", "Chitin", "Glycogen", "Lipid", "Protein", "Dryweight"))) %>%
arrange(vars, Parameter) %>%
select(Parameter, Estimate = est, `Est. error` = err,
`CI lower` = lwr, `CI upper` = upr, `p` = p_val, star) %>%
rename(` ` = star) %>%
#write_csv('output/metabolite_SEM_table.csv')
mutate(p = ifelse(p > 0.001, round(p, 3), '< 0.001')) %>%
kable() %>%
kable_styling(full_width = FALSE) %>%
group_rows("Carbohydrates", 1, 4) %>%
group_rows("Chitin", 5, 8) %>%
group_rows("Glycogen", 9, 12) %>%
group_rows("Lipids", 13, 16) %>%
group_rows("Protein", 17, 20) %>%
group_rows("Dry weight", 21, 23)| Parameter | Estimate | Est. error | CI lower | CI upper | p | |
|---|---|---|---|---|---|---|
| Carbohydrates | ||||||
| Dry weight | 0.087 | 0.265 | -0.441 | 0.600 | 0.364 | |
| Sex (M) | 0.004 | 0.419 | -0.804 | 0.844 | 0.499 | |
| Treatment (E) | -0.245 | 0.303 | -0.838 | 0.353 | 0.205 | |
| Sex (M) x Treatment (E) | -0.418 | 0.348 | -1.101 | 0.253 | 0.117 | |
| Chitin | ||||||
| Dry weight | -0.484 | 0.268 | -1.014 | 0.041 | 0.036 | * |
| Sex (M) | 0.398 | 0.430 | -0.455 | 1.244 | 0.176 | |
| Treatment (E) | -0.120 | 0.282 | -0.658 | 0.446 | 0.326 | |
| Sex (M) x Treatment (E) | -0.313 | 0.325 | -0.962 | 0.328 | 0.161 | |
| Glycogen | ||||||
| Dry weight | 0.339 | 0.263 | -0.177 | 0.856 | 0.101 | |
| Sex (M) | -0.268 | 0.419 | -1.094 | 0.550 | 0.259 | |
| Treatment (E) | 0.216 | 0.297 | -0.359 | 0.803 | 0.233 | |
| Sex (M) x Treatment (E) | 0.394 | 0.346 | -0.292 | 1.059 | 0.126 | |
| Lipids | ||||||
| Dry weight | 0.537 | 0.255 | 0.026 | 1.028 | 0.019 | * |
| Sex (M) | -0.120 | 0.415 | -0.930 | 0.680 | 0.393 | |
| Treatment (E) | 0.392 | 0.271 | -0.152 | 0.918 | 0.075 | |
| Sex (M) x Treatment (E) | -0.053 | 0.314 | -0.681 | 0.547 | 0.431 | |
| Protein | ||||||
| Dry weight | -0.211 | 0.271 | -0.734 | 0.331 | 0.214 | |
| Sex (M) | -0.009 | 0.427 | -0.831 | 0.828 | 0.491 | |
| Treatment (E) | -0.211 | 0.312 | -0.826 | 0.402 | 0.251 | |
| Sex (M) x Treatment (E) | 0.343 | 0.357 | -0.366 | 1.045 | 0.166 | |
| Dry weight | ||||||
| Sex (M) | -1.616 | 0.140 | -1.894 | -1.347 | < 0.001 | * |
| Treatment (E) | 0.522 | 0.156 | 0.211 | 0.821 | 0.002 | * |
| Sex (M) x Treatment (E) | -0.352 | 0.196 | -0.741 | 0.043 | 0.038 | * |
Complete output from summary.brmsfit()
- ‘Group-Level Effects’ (also called random effects): This shows the (co)variances associated with the line-specific intercepts (which have names like
sd(Lipid_Intercept)) and slopes (e.g.sd(Dryweight_treatmentPolyandry)), as well as the correlations between these effects (e.g.cor(Lipid_Intercept,Protein_Intercept)is the correlation in line effects on lipids and proteins) - ‘Population-Level Effects:’ (also called fixed effects): These give the estimates of the intercept (i.e. for female M flies) and the effects of treatment, sex, dry weight, and the treatment \(\times\) sex interaction, for each response variable.
- ‘Family Specific Parameters’: This is the parameter sigma for the residual variance for each response variable
- ‘Residual Correlations:’ This give the correlations between the residuals for each pairs of response variables.
Note that the model has converged (Rhat = 1) and the posterior is adequately samples (high ESS values).
brms_metabolite_SEM Family: MV(gaussian, gaussian, gaussian, gaussian, gaussian, gaussian)
Links: mu = identity; sigma = identity
mu = identity; sigma = identity
mu = identity; sigma = identity
mu = identity; sigma = identity
mu = identity; sigma = identity
mu = identity; sigma = identity
Formula: Lipid ~ sex * treatment + Dryweight + (1 | p | line)
Carbohydrate ~ sex * treatment + Dryweight + (1 | p | line)
Protein ~ sex * treatment + Dryweight + (1 | p | line)
Glycogen ~ sex * treatment + Dryweight + (1 | p | line)
Chitin ~ sex * treatment + Dryweight + (1 | p | line)
Dryweight ~ sex * treatment + (1 | p | line)
Data: scaled_metabolites %>% rename(Dryweight = Dry_weig (Number of observations: 48)
Draws: 4 chains, each with iter = 5000; warmup = 2500; thin = 1;
total post-warmup draws = 10000
Group-Level Effects:
~line (Number of levels: 8)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Lipid_Intercept) 0.09 0.15 0.00 0.50 1.00 2226 3040
sd(Carbohydrate_Intercept) 0.07 0.13 0.00 0.50 1.00 3422 2717
sd(Protein_Intercept) 0.05 0.12 0.00 0.44 1.00 4846 3055
sd(Glycogen_Intercept) 0.03 0.07 0.00 0.26 1.00 9019 5018
sd(Chitin_Intercept) 0.07 0.12 0.00 0.44 1.00 3130 2735
sd(Dryweight_Intercept) 0.05 0.08 0.00 0.27 1.00 2860 4008
cor(Lipid_Intercept,Carbohydrate_Intercept) 0.02 0.38 -0.70 0.71 1.00 13199 6832
cor(Lipid_Intercept,Protein_Intercept) 0.00 0.38 -0.70 0.71 1.00 14787 6474
cor(Carbohydrate_Intercept,Protein_Intercept) 0.01 0.38 -0.70 0.72 1.00 11274 6759
cor(Lipid_Intercept,Glycogen_Intercept) 0.01 0.37 -0.69 0.71 1.00 15679 6795
cor(Carbohydrate_Intercept,Glycogen_Intercept) 0.01 0.37 -0.69 0.70 1.00 11674 7312
cor(Protein_Intercept,Glycogen_Intercept) 0.00 0.38 -0.71 0.71 1.00 9312 7759
cor(Lipid_Intercept,Chitin_Intercept) 0.00 0.38 -0.70 0.71 1.00 12344 7086
cor(Carbohydrate_Intercept,Chitin_Intercept) -0.00 0.38 -0.71 0.71 1.00 11084 7659
cor(Protein_Intercept,Chitin_Intercept) 0.01 0.38 -0.70 0.71 1.00 8263 7879
cor(Glycogen_Intercept,Chitin_Intercept) 0.01 0.38 -0.70 0.71 1.00 7362 7966
cor(Lipid_Intercept,Dryweight_Intercept) -0.01 0.37 -0.69 0.69 1.00 11583 7502
cor(Carbohydrate_Intercept,Dryweight_Intercept) 0.01 0.38 -0.70 0.72 1.00 9377 7379
cor(Protein_Intercept,Dryweight_Intercept) 0.01 0.38 -0.70 0.71 1.00 8299 7795
cor(Glycogen_Intercept,Dryweight_Intercept) -0.01 0.38 -0.70 0.71 1.00 7296 8107
cor(Chitin_Intercept,Dryweight_Intercept) -0.01 0.38 -0.71 0.70 1.00 7056 8688
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Lipid_Intercept -0.12 0.24 -0.58 0.36 1.00 8283 7353
Carbohydrate_Intercept 0.22 0.27 -0.31 0.77 1.00 9902 7985
Protein_Intercept 0.02 0.28 -0.52 0.56 1.00 10823 8268
Glycogen_Intercept -0.07 0.26 -0.58 0.43 1.00 10745 8442
Chitin_Intercept -0.06 0.26 -0.57 0.45 1.00 8522 6813
Dryweight_Intercept 0.63 0.11 0.42 0.85 1.00 8348 7726
Lipid_sexMale -0.12 0.42 -0.93 0.68 1.00 6798 7196
Lipid_treatmentPolyandry 0.39 0.27 -0.15 0.92 1.00 8798 7658
Lipid_Dryweight 0.54 0.25 0.03 1.03 1.00 5568 6453
Lipid_sexMale:treatmentPolyandry -0.05 0.31 -0.68 0.55 1.00 9816 7636
Carbohydrate_sexMale 0.00 0.42 -0.80 0.84 1.00 7962 7496
Carbohydrate_treatmentPolyandry -0.24 0.30 -0.84 0.35 1.00 10245 7539
Carbohydrate_Dryweight 0.09 0.26 -0.44 0.60 1.00 6631 7408
Carbohydrate_sexMale:treatmentPolyandry -0.42 0.35 -1.10 0.25 1.00 11374 7608
Protein_sexMale -0.01 0.43 -0.83 0.83 1.00 8339 7487
Protein_treatmentPolyandry -0.21 0.31 -0.83 0.40 1.00 10374 8676
Protein_Dryweight -0.21 0.27 -0.73 0.33 1.00 6785 7231
Protein_sexMale:treatmentPolyandry 0.34 0.36 -0.37 1.04 1.00 10331 7256
Glycogen_sexMale -0.27 0.42 -1.09 0.55 1.00 7302 7487
Glycogen_treatmentPolyandry 0.22 0.30 -0.36 0.80 1.00 9860 8345
Glycogen_Dryweight 0.34 0.26 -0.18 0.86 1.00 6689 7518
Glycogen_sexMale:treatmentPolyandry 0.39 0.35 -0.29 1.06 1.00 11426 8369
Chitin_sexMale 0.40 0.43 -0.45 1.24 1.00 7576 7166
Chitin_treatmentPolyandry -0.12 0.28 -0.66 0.45 1.00 8551 7564
Chitin_Dryweight -0.48 0.27 -1.01 0.04 1.00 6202 6889
Chitin_sexMale:treatmentPolyandry -0.31 0.32 -0.96 0.33 1.00 9734 8005
Dryweight_sexMale -1.62 0.14 -1.89 -1.35 1.00 8443 7594
Dryweight_treatmentPolyandry 0.52 0.16 0.21 0.82 1.00 6327 7005
Dryweight_sexMale:treatmentPolyandry -0.35 0.20 -0.74 0.04 1.00 7145 7157
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma_Lipid 0.74 0.09 0.59 0.93 1.00 6746 7390
sigma_Carbohydrate 1.00 0.11 0.81 1.25 1.00 9278 7809
sigma_Protein 1.02 0.12 0.82 1.28 1.00 9185 7780
sigma_Glycogen 0.95 0.11 0.76 1.19 1.00 11909 7770
sigma_Chitin 0.85 0.10 0.67 1.06 1.00 7486 7023
sigma_Dryweight 0.36 0.04 0.29 0.46 1.00 8239 7406
Residual Correlations:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
rescor(Lipid,Carbohydrate) -0.32 0.14 -0.58 -0.04 1.00 5457 6371
rescor(Lipid,Protein) -0.06 0.15 -0.34 0.24 1.00 10423 7214
rescor(Carbohydrate,Protein) 0.04 0.14 -0.25 0.31 1.00 11539 7956
rescor(Lipid,Glycogen) 0.04 0.15 -0.26 0.33 1.00 7796 6676
rescor(Carbohydrate,Glycogen) -0.01 0.15 -0.29 0.28 1.00 11161 7736
rescor(Protein,Glycogen) -0.13 0.14 -0.41 0.16 1.00 10372 7470
rescor(Lipid,Chitin) -0.07 0.15 -0.35 0.22 1.00 8329 7466
rescor(Carbohydrate,Chitin) -0.42 0.12 -0.64 -0.16 1.00 9965 8172
rescor(Protein,Chitin) 0.06 0.15 -0.22 0.35 1.00 9097 7373
rescor(Glycogen,Chitin) -0.04 0.15 -0.32 0.25 1.00 9577 7869
rescor(Lipid,Dryweight) 0.16 0.18 -0.20 0.49 1.00 5978 6731
rescor(Carbohydrate,Dryweight) 0.02 0.17 -0.33 0.35 1.00 5275 5847
rescor(Protein,Dryweight) -0.04 0.17 -0.38 0.30 1.00 6616 7322
rescor(Glycogen,Dryweight) 0.01 0.17 -0.33 0.35 1.00 6280 7242
rescor(Chitin,Dryweight) 0.25 0.18 -0.12 0.57 1.00 4907 5998
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
Posterior effect size of treatment on metabolite abundance, for each sex
Here, we use the model to predict the mean concentration of each metabolite (in standard units) in each treatment and sex (averaged across the eight replicate selection lines). We then calculate the effect size of treatment by subtracting the (sex-specific) mean for the M treatment from the mean for the E treatment; thus a value of 1 would mean that the E treatment has a mean that is larger by 1 standard deviation. Thus, the y-axis in the following graphs essentially shows the posterior estimate of standardised effect size (Cohen’s d), from the model shown above.
Because the model contains dry weight as a mediator variable, we created these predictions two different ways, and display the answer for both using tabs in the following figures/tables. Firstly, we predicted the means controlling for differences in dry weight between sexes and treatments; this was done by deriving the predictions with dry weight set to its global mean, for both sexes and treatments. Secondly, we derived predictions without controlling for dry weight. This was done by deriving the predictions with dry weight set to its average value for the appropriate treatment-sex combination.
By clicking the tabs and comparing, one can see that the estimates of the treatment effect hardly change when differences in dry weight are controlled for. This indicates that dry mass does not have an important role in mediating the effect of treatment on metabolite composition, even though body size differs between treatments. Thus, we conclude that the M vs E treatments caused metabolite composition to evolve, through mechanisms other than the evolution of dry weight.
Figure
Not controlling for differences in dry weight between treatments
new <- expand_grid(sex = c("Male", "Female"),
treatment = c("Monogamy", "Polyandry"),
Dryweight = NA, line = NA) %>%
mutate(type = 1:n())
levels <- c("Carbohydrate", "Chitin", "Glycogen", "Lipid", "Protein", "Dryweight")
# Estimate mean dry weight for each of the 4 sex/treatment combinations
evolved_mean_dryweights <- data.frame(
new[,1:2],
fitted(brms_metabolite_SEM, re_formula = NA,
newdata = new %>% select(-Dryweight),
summary = TRUE, resp = "Dryweight")) %>%
as_tibble()
# Find the mean dry weight for males and females (across treatments)
male_dryweight <- mean(evolved_mean_dryweights$Estimate[1:2])
female_dryweight <- mean(evolved_mean_dryweights$Estimate[3:4])
new_metabolites <- bind_rows(
expand_grid(sex = c("Male", "Female"),
treatment = c("Monogamy", "Polyandry"),
Dryweight = c(male_dryweight, female_dryweight), line = NA) %>%
filter(sex == "Male" & Dryweight == male_dryweight |
sex == "Female" & Dryweight == female_dryweight) %>%
mutate(type = 1:4),
evolved_mean_dryweights %>% select(sex, treatment, Dryweight = Estimate) %>%
mutate(line = NA, type = 5:8)
)
# Predict data from the SEM of metabolites...
# Because we use sum contrasts for "line" and line=NA in the new data,
# this function predicts at the global means across the 4 lines (see ?posterior_epred)
fitted_values <- posterior_epred(
brms_metabolite_SEM, newdata = new_metabolites, re_formula = NA,
summary = FALSE) %>%
reshape2::melt() %>% rename(draw = Var1, type = Var2, variable = Var3) %>%
as_tibble() %>%
left_join(new_metabolites, by = "type") %>%
select(draw, variable, value, sex, treatment, Dryweight) %>%
mutate(variable = factor(variable, levels))
treat_diff_standard_dryweight <- fitted_values %>%
filter(Dryweight %in% c(male_dryweight, female_dryweight)) %>%
spread(treatment, value) %>%
mutate(`Difference in means (Poly - Mono)` = Polyandry - Monogamy)
treat_diff_actual_dryweight <- fitted_values %>%
filter(!(Dryweight %in% c(male_dryweight, female_dryweight))) %>%
select(-Dryweight) %>%
spread(treatment, value) %>%
mutate(`Difference in means (Poly - Mono)` = Polyandry - Monogamy)
summary_dat1 <- treat_diff_actual_dryweight %>%
filter(variable != 'Dryweight') %>%
rename(x = `Difference in means (Poly - Mono)`) %>%
group_by(variable, sex) %>%
summarise(`Difference in means (Poly - Mono)` = median(x),
`Lower 95% CI` = quantile(x, probs = 0.025),
`Upper 95% CI` = quantile(x, probs = 0.975),
p = 1 - as.numeric(bayestestR::p_direction(x)),
` ` = ifelse(p < 0.05, "\\*", ""),
.groups = "drop")
summary_dat2 <- treat_diff_standard_dryweight %>%
filter(variable != 'Dryweight') %>%
rename(x = `Difference in means (Poly - Mono)`) %>%
group_by(variable, sex) %>%
summarise(`Difference in means (Poly - Mono)` = median(x),
`Lower 95% CI` = quantile(x, probs = 0.025),
`Upper 95% CI` = quantile(x, probs = 0.975),
p = 1 - as.numeric(bayestestR::p_direction(x)),
` ` = ifelse(p < 0.05, "\\*", ""),
.groups = "drop")
sampled_draws <- sample(unique(fitted_values$draw), 100)
ylims <- c(-1.8, 1.8)
p1 <- treat_diff_actual_dryweight %>%
filter(variable != 'Dryweight') %>%
ggplot(aes(x = sex, y = `Difference in means (Poly - Mono)`,fill = sex)) +
geom_hline(yintercept = 0, linetype = 2) +
stat_halfeye() +
# geom_line(data = treat_diff_actual_dryweight %>%
# filter(draw %in% sampled_draws) %>%
# filter(variable != 'Dryweight'),
# alpha = 0.8, size = 0.12, colour = "black", aes(group = draw)) +
geom_point(data = summary_dat1, pch = 21, colour = "black", size = 3.1) +
scale_fill_brewer(palette = 'Pastel1', direction = 1, name = "") +
scale_colour_brewer(palette = 'Pastel1', direction = 1, name = "") +
facet_wrap( ~ variable, nrow = 1) +
theme_bw() +
theme(legend.position = 'none',
strip.background = element_blank(),
text = element_text(family = nice_font),
panel.grid.major.x = element_blank()) +
coord_cartesian(ylim = ylims) +
ylab("Difference in means between\nselection treatments (E - M)") + xlab("Sex") +
#ggsave("figures/metabolite_plot_nolines.pdf", height=2.5, width=6) +
NULL
p1
| Version | Author | Date |
|---|---|---|
| aa3e11e | MartinGarlovsky | 2021-10-05 |
| 21567bb | MartinGarlovsky | 2021-03-12 |
| ffb09dd | MartinGarlovsky | 2021-02-08 |
| 034431b | lukeholman | 2020-12-18 |
| b6cf554 | lukeholman | 2020-12-11 |
| a88f037 | lukeholman | 2020-12-11 |
| bb96acf | lukeholman | 2020-12-11 |
| e5c580f | lukeholman | 2020-12-11 |
| 7fca240 | lukeholman | 2020-12-10 |
| f7c88a2 | lukeholman | 2020-12-10 |
| 43cc270 | lukeholman | 2020-12-09 |
Figure XX: Posterior estimates of the treatment effect size for both sexes, for each of the five metabolites. A positive value means that the mean metabolite concentration is higher in the E treatment than the M treatment, while a negative effects denotes M > E. A strongly supported treatment effect is implied by the majority of the posterior lying to one side of zero. The error bars summarise the 66% and 95% quantiles of the posterior. This plot was created used posterior predictions of the means that were not adjusted for differences in dry weight between treatments.
Controlling for differences in dry weight between treatments
treat_diff_standard_dryweight %>%
filter(variable != 'Dryweight') %>%
ggplot(aes(x = sex, y = `Difference in means (Poly - Mono)`, fill = sex)) +
geom_hline(yintercept = 0, linetype = 2) +
stat_halfeye() +
# geom_line(data = treat_diff_standard_dryweight %>%
# filter(draw %in% sampled_draws) %>%
# filter(variable != 'Dryweight'),
# alpha = 0.8, size = 0.12, colour = "black", aes(group = draw)) +
geom_point(data = summary_dat2, pch = 21, colour = "black", size = 3.1) +
scale_fill_brewer(palette = 'Pastel1', direction = 1, name = "") +
scale_colour_brewer(palette = 'Pastel1', direction = 1, name = "") +
facet_wrap( ~ variable, nrow = 1) +
theme_bw() +
theme(legend.position = 'none',
strip.background = element_blank(),
text = element_text(family = nice_font),
panel.grid.major.x = element_blank()) +
coord_cartesian(ylim = ylims) +
ylab("Difference in means between\nselection treatments (E - M)") + xlab("Sex") +
#ggsave('figures/metabolite_plotCONTROLLED_nolines.pdf', height=2.5, width=6) +
NULL
| Version | Author | Date |
|---|---|---|
| aa3e11e | MartinGarlovsky | 2021-10-05 |
| 21567bb | MartinGarlovsky | 2021-03-12 |
| ffb09dd | MartinGarlovsky | 2021-02-08 |
| 034431b | lukeholman | 2020-12-18 |
| b6cf554 | lukeholman | 2020-12-11 |
| a88f037 | lukeholman | 2020-12-11 |
| bb96acf | lukeholman | 2020-12-11 |
| e5c580f | lukeholman | 2020-12-11 |
| 7fca240 | lukeholman | 2020-12-10 |
| f7c88a2 | lukeholman | 2020-12-10 |
Figure XX: Posterior estimates of the treatment effect size for both sexes, for each of the five metabolites. A positive value means that the mean metabolite concentration is higher in the P treatment than the M treatment, while a negative effects denotes M > E. A strongly supported treatment effect is implied by the majority of the posterior lying to one side of zero. The error bars summarise the 66% and 95% quantiles of the posterior. This plot was created used posterior predictions of the means that were adjusted for differences in dry weight between treatments.
Table
Not controlling for differences in dry weight between treatments
summary_dat1 %>%
kable(digits = 3) %>%
kable_styling(full_width = FALSE)| variable | sex | Difference in means (Poly - Mono) | Lower 95% CI | Upper 95% CI | p | |
|---|---|---|---|---|---|---|
| Carbohydrate | Female | -0.198 | -0.759 | 0.370 | 0.237 | |
| Carbohydrate | Male | -0.650 | -1.320 | 0.034 | 0.030 | * |
| Chitin | Female | -0.374 | -0.884 | 0.143 | 0.077 | |
| Chitin | Male | -0.518 | -1.087 | 0.065 | 0.041 | * |
| Glycogen | Female | 0.391 | -0.139 | 0.931 | 0.076 | |
| Glycogen | Male | 0.672 | 0.035 | 1.285 | 0.021 | * |
| Lipid | Female | 0.674 | 0.184 | 1.150 | 0.006 | * |
| Lipid | Male | 0.432 | -0.139 | 0.979 | 0.064 | |
| Protein | Female | -0.321 | -0.899 | 0.268 | 0.138 | |
| Protein | Male | 0.101 | -0.575 | 0.773 | 0.386 |
Controlling for differences in dry weight between treatments
summary_dat2 %>%
kable(digits = 3) %>%
kable_styling(full_width = FALSE)| variable | sex | Difference in means (Poly - Mono) | Lower 95% CI | Upper 95% CI | p | |
|---|---|---|---|---|---|---|
| Carbohydrate | Female | -0.245 | -0.838 | 0.353 | 0.205 | |
| Carbohydrate | Male | -0.663 | -1.326 | 0.015 | 0.027 | * |
| Chitin | Female | -0.124 | -0.658 | 0.446 | 0.326 | |
| Chitin | Male | -0.435 | -1.002 | 0.147 | 0.074 | |
| Glycogen | Female | 0.214 | -0.359 | 0.803 | 0.233 | |
| Glycogen | Male | 0.612 | -0.023 | 1.225 | 0.029 | * |
| Lipid | Female | 0.390 | -0.152 | 0.918 | 0.075 | |
| Lipid | Male | 0.343 | -0.217 | 0.888 | 0.113 | |
| Protein | Female | -0.210 | -0.826 | 0.402 | 0.251 | |
| Protein | Male | 0.132 | -0.539 | 0.803 | 0.343 |
Posterior difference in treatment effect size between sexes
This section essentially examines the treatment \(\times\) sex interaction term, by calculating the difference in the effect size of the E/M treatment between sexes, for each of the five metabolites. We find no strong evidence for a treatment \(\times\) sex interaction, i.e. the treatment effects did not differ detectably between sexes.
Figure
Not controlling for differences in dry weight between treatments
treatsex_interaction_data1 <- treat_diff_actual_dryweight %>%
select(draw, variable, sex, d = `Difference in means (Poly - Mono)`) %>%
arrange(draw, variable, sex) %>%
group_by(draw, variable) %>%
summarise(`Difference in effect size between sexes (male - female)` = d[2] - d[1],
.groups = "drop") # males - females
p2 <- treatsex_interaction_data1 %>%
filter(variable != 'Dryweight') %>%
ggplot(aes(x = `Difference in effect size between sexes (male - female)`, y = 1, fill = stat(x < 0))) +
geom_vline(xintercept = 0, linetype = 2) +
stat_halfeye() +
facet_wrap( ~ variable) +
scale_fill_brewer(palette = 'Pastel2', direction = 1, name = "") +
theme_bw() +
theme(legend.position = 'none',
text = element_text(family = nice_font),
strip.background = element_blank()) +
ylab("Posterior density") +
#ggsave("figures/metabolite_interaction_plot.pdf", height=4, width=6) +
NULL
p2
Figure XX: Posterior estimates of the difference in the treatment effect size (i.e. mean of E minus mean of M) between males and females, for each of the five metabolites. A positive value means that the effect size is more positive in males, and negative means it is more positive in females. A strongly supported sex difference in effect size would be implied by the majority of the posterior lying to one side of zero. The error bars summarise the 66% and 95% quantiles of the posterior. This plot was created used posterior predictions of the means that were not adjusted for differences in dry weight between treatments.
Controlling for differences in dry weight between treatments
treatsex_interaction_data2 <- treat_diff_standard_dryweight %>%
select(draw, variable, sex, d = `Difference in means (Poly - Mono)`) %>%
arrange(draw, variable, sex) %>%
group_by(draw, variable) %>%
summarise(`Difference in effect size between sexes (male - female)` = d[2] - d[1],
.groups = "drop") # males - females
treatsex_interaction_data2 %>%
filter(variable != 'Dryweight') %>%
ggplot(aes(x = `Difference in effect size between sexes (male - female)`, y = 1, fill = stat(x < 0))) +
geom_vline(xintercept = 0, linetype = 2) +
stat_halfeye() +
facet_wrap( ~ variable) +
scale_fill_brewer(palette = 'Pastel2', direction = 1, name = "") +
theme_bw() +
theme(legend.position = 'none',
text = element_text(family = nice_font),
strip.background = element_blank()) +
ylab("Posterior density")
Figure XX: Posterior estimates of the difference in the treatment effect size (i.e. mean of E minus mean of M) between males and females, for each of the five metabolites. A positive value means that the effect size is more positive in males, and negative means it is more positive in females. A strongly supported sex difference in effect size would be implied by the majority of the posterior lying to one side of zero. The error bars summarise the 66% and 95% quantiles of the posterior. This plot was created used posterior predictions of the means that were adjusted for differences in dry weight between treatments.
Table
Not controlling for differences in dry weight between treatments
treatsex_interaction_data1 %>%
filter(variable != 'Dryweight') %>%
rename(x = `Difference in effect size between sexes (male - female)`) %>%
group_by(variable) %>%
summarise(`Difference in effect size between sexes (male - female)` = median(x),
`Lower 95% CI` = quantile(x, probs = 0.025),
`Upper 95% CI` = quantile(x, probs = 0.975),
p = 1 - as.numeric(bayestestR::p_direction(x)),
` ` = ifelse(p < 0.05, "\\*", ""),
.groups = "drop") %>%
kable(digits=3) %>%
kable_styling(full_width = FALSE)| variable | Difference in effect size between sexes (male - female) | Lower 95% CI | Upper 95% CI | p | |
|---|---|---|---|---|---|
| Carbohydrate | -0.451 | -1.103 | 0.209 | 0.092 | |
| Chitin | -0.139 | -0.743 | 0.479 | 0.318 | |
| Glycogen | 0.279 | -0.381 | 0.900 | 0.202 | |
| Lipid | -0.244 | -0.831 | 0.326 | 0.207 | |
| Protein | 0.422 | -0.255 | 1.098 | 0.112 |
Controlling for differences in dry weight between treatments
treatsex_interaction_data2 %>%
filter(variable != 'Dryweight') %>%
rename(x = `Difference in effect size between sexes (male - female)`) %>%
group_by(variable) %>%
summarise(`Difference in effect size between sexes (male - female)` = median(x),
`Lower 95% CI` = quantile(x, probs = 0.025),
`Upper 95% CI` = quantile(x, probs = 0.975),
p = 1 - as.numeric(bayestestR::p_direction(x)),
` ` = ifelse(p < 0.05, "\\*", ""),
.groups = "drop") %>%
kable(digits=3) %>%
kable_styling(full_width = FALSE)| variable | Difference in effect size between sexes (male - female) | Lower 95% CI | Upper 95% CI | p | |
|---|---|---|---|---|---|
| Carbohydrate | -0.420 | -1.101 | 0.253 | 0.117 | |
| Chitin | -0.308 | -0.962 | 0.328 | 0.161 | |
| Glycogen | 0.399 | -0.292 | 1.059 | 0.126 | |
| Lipid | -0.053 | -0.681 | 0.547 | 0.431 | |
| Protein | 0.344 | -0.366 | 1.045 | 0.166 |
sessionInfo()R version 4.0.3 (2020-10-10) Platform: x86_64-apple-darwin17.0 (64-bit) Running under: macOS Big Sur 10.16 Matrix products: default BLAS: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib locale: [1] en_GB.UTF-8/en_GB.UTF-8/en_GB.UTF-8/C/en_GB.UTF-8/en_GB.UTF-8 attached base packages: [1] stats graphics grDevices utils datasets methods base other attached packages: [1] showtext_0.9-4 showtextdb_3.0 sysfonts_0.8.5 knitrhooks_0.0.4 knitr_1.33 kableExtra_1.3.4 DT_0.18 tidybayes_3.0.1 brms_2.16.1 [10] Rcpp_1.0.7 ggridges_0.5.3 gridExtra_2.3 GGally_2.1.2 forcats_0.5.1 stringr_1.4.0 dplyr_1.0.7 purrr_0.3.4 readr_2.0.1 [19] tidyr_1.1.3 tibble_3.1.3 ggplot2_3.3.5 tidyverse_1.3.1 workflowr_1.6.2 loaded via a namespace (and not attached): [1] readxl_1.3.1 backports_1.2.1 systemfonts_1.0.2 plyr_1.8.6 igraph_1.2.6 svUnit_1.0.6 splines_4.0.3 [8] crosstalk_1.1.1 rstantools_2.1.1 inline_0.3.19 digest_0.6.27 htmltools_0.5.1.1 rsconnect_0.8.24 fansi_0.5.0 [15] magrittr_2.0.1 checkmate_2.0.0 tzdb_0.1.2 modelr_0.1.8 RcppParallel_5.1.4 matrixStats_0.60.0 vroom_1.5.4 [22] svglite_2.0.0 xts_0.12.1 prettyunits_1.1.1 colorspace_2.0-2 rvest_1.0.1 ggdist_3.0.0 haven_2.4.3 [29] xfun_0.25 callr_3.7.0 crayon_1.4.1 jsonlite_1.7.2 lme4_1.1-27.1 zoo_1.8-9 glue_1.4.2 [36] gtable_0.3.0 webshot_0.5.2 V8_3.4.2 distributional_0.2.2 pkgbuild_1.2.0 rstan_2.21.1 abind_1.4-5 [43] scales_1.1.1 mvtnorm_1.1-2 DBI_1.1.1 miniUI_0.1.1.1 viridisLite_0.4.0 xtable_1.8-4 bit_4.0.4 [50] stats4_4.0.3 StanHeaders_2.21.0-7 datawizard_0.2.0 htmlwidgets_1.5.3 httr_1.4.2 DiagrammeR_1.0.6.1 threejs_0.3.3 [57] arrayhelpers_1.1-0 RColorBrewer_1.1-2 posterior_1.0.1 ellipsis_0.3.2 pkgconfig_2.0.3 reshape_0.8.8 loo_2.4.1 [64] farver_2.1.0 sass_0.4.0 dbplyr_2.1.1 utf8_1.2.2 labeling_0.4.2 tidyselect_1.1.1 rlang_0.4.11 [71] reshape2_1.4.4 later_1.3.0 visNetwork_2.0.9 munsell_0.5.0 cellranger_1.1.0 tools_4.0.3 cli_3.0.1 [78] generics_0.1.0 broom_0.7.9 evaluate_0.14 fastmap_1.1.0 yaml_2.2.1 bit64_4.0.5 processx_3.5.2 [85] fs_1.5.0 nlme_3.1-152 whisker_0.4 mime_0.11 projpred_2.0.2 xml2_1.3.2 compiler_4.0.3 [92] bayesplot_1.8.1 shinythemes_1.2.0 rstudioapi_0.13 gamm4_0.2-6 curl_4.3.2 reprex_2.0.1 bslib_0.2.5.1 [99] stringi_1.7.3 highr_0.9 ps_1.6.0 Brobdingnag_1.2-6 lattice_0.20-44 Matrix_1.3-4 nloptr_1.2.2.2 [106] markdown_1.1 shinyjs_2.0.0 tensorA_0.36.2 vctrs_0.3.8 pillar_1.6.2 lifecycle_1.0.0 jquerylib_0.1.4 [113] bridgesampling_1.1-2 insight_0.14.3 httpuv_1.6.2 R6_2.5.1 promises_1.2.0.1 codetools_0.2-18 boot_1.3-28 [120] colourpicker_1.1.0 MASS_7.3-54 gtools_3.9.2 assertthat_0.2.1 rprojroot_2.0.2 withr_2.4.2 shinystan_2.5.0 [127] bayestestR_0.10.5 mgcv_1.8-36 parallel_4.0.3 hms_1.1.0 grid_4.0.3 coda_0.19-4 minqa_1.2.4 [134] rmarkdown_2.10 git2r_0.28.0 shiny_1.6.0 lubridate_1.7.10 base64enc_0.1-3 dygraphs_1.1.1.6