knitr::opts_chunk$set(echo = TRUE) knitr::opts_chunk$set(message = FALSE) knitr::opts_chunk$set(warning = FALSE) library(ProvidenciaChemo) library(tidyverse) library(patchwork) library(tidybayes) library(modelr) library(emmeans) library(magrittr) library(emmeans) library(rstanarm) theme_set(theme_my)
tyrosine_data <- read_csv(here::here('extdata','Figure_Ext7a_LTyr.csv')) %>% mutate(food = fct_relevel(food, "OP50"), cond = fct_relevel(cond, "off_food"), data_type = "raw") %>% filter(date != "2018_10_23") #too few JUb39 measured to include these dates stan_glm <- rstanarm::stan_lmer(tyrosine_data, formula = log(response.time) ~ food*cond + (food | date) + (1|plateID), seed = 6994, cores = 4, chains = 4) #-----------posterior HDI estimates for relative TA effects--------- tyrosine_data <- tyrosine_data %>% # mutate(response.time = case_when( # response.time < 1 ~ 1, # TRUE ~ response.time)) %>% group_by(cond, food, date) %>% summarise(meanTA_0 = mean(log(response.time))) %>% filter(food == "OP50") %>% ungroup() %>% select(date, cond, meanTA_0) %>% full_join(., tyrosine_data) %>% mutate(rel_log = log(response.time) - meanTA_0, food = fct_relevel(food, "OP50"), cond = fct_relevel(cond, "off_food")) contrasts.95 <- modelbased::estimate_contrasts(stan_glm, ci = 0.95)[c(1,6),] %>% mutate(food = factor("JUb39", levels = c("OP50", "JUb39")), cond = c("off_food", "Tyrosine")) %>% select(food, cond, m = Median, ll = CI_low, hh = CI_high) contrasts.66 <- modelbased::estimate_contrasts(stan_glm, ci = 0.66)[c(1,6),] %>% mutate(food = factor("JUb39", levels = c("OP50", "JUb39")), cond = c("off_food", "Tyrosine")) %>% select(food, cond, m = Median, l = CI_low, h = CI_high) fitted <- full_join(contrasts.95, contrasts.66) %>% mutate(food = fct_relevel(food, "OP50"), cond = fct_relevel(cond, "off_food")) %>% mutate_if(is.numeric, function(.) {-1 * .}) # emmeans(stan_glm, pairwise ~ (cond| genotype))$contrasts %>% # coda::as.mcmc() %>% # bayesplot::mcmc_intervals_data(prob = 0.66, prob_outer = 0.95) %>% # mutate(genotype = c("N2", "octr-1"), # data_type = "fit", # cond = "Tyr+10mM_TA") %>% # mutate(cond = factor(cond, levels = c("Tyrosine", "Tyr+10mM_TA")), # genotype = fct_relevel(genotype, c("N2", "octr-1"))) %>% # mutate_if(is.numeric, funs(. * -1)) #plot plot <- tyrosine_data %>% ggplot(aes(x = data_type)) + geom_relLatency(fitted = fitted, fillvar = food, dotvar = food, yvar = rel_log) + scale_color_plot("grey-blue", drop = TRUE) + scale_fill_plot("grey-blue", drop = TRUE) + facet_grid(.~cond+food) + add.n('data_type', y.pos = -1.6) + labs(x = "genotype", y = "relative reversal latency [log(s)]") + guides(colour = FALSE, fill = FALSE) + theme(axis.text.x = element_blank()) plot #------------posterior draws for non-normalized data------- fitted <- tyrosine_data %>% data_grid(food,cond) %>% add_fitted_draws(stan_glm, re_formula = NA) %>% mutate(response.time = exp(.value), data_type = "fit") plot <- tyrosine_data %>% mutate(response.time = case_when( response.time < 1 ~ 1, TRUE ~ response.time )) %>% ggplot(aes(x = data_type, y = response.time)) + stat_summary(geom = "bar", fun.y = mean, aes(fill = food), width = 0.6, alpha = 0.5) + ggbeeswarm::geom_quasirandom(aes(colour = food), width = 0.2, alpha = 0.75) + facet_grid(.~cond + food, labeller = labeller(genotype = label_wrap_gen(10), food = as_labeller(c("OP50" = "", "JUb39" = "")))) + scale_color_plot(palette = "grey-blue", drop = TRUE) + scale_fill_plot(palette = "grey-blue", drop = TRUE) + stat_summary(geom = "errorbar", fun.data = mean_se, width = 0.2) + add.n('data_type', y.pos = 0.9) + stat_pointinterval(aes(y=response.time, x = 1.5), data = fitted, fatten_point = 0, size_range = c(0.3, 1), colour = "darkgrey") + stat_summary(data = fitted, aes(y=response.time, x = 1.5), fun.y = median, fun.ymin = median, fun.ymax = median, geom = "crossbar", width = 0.05, lwd = 0.35, colour = "darkgrey") + labs(x = "genotype", y = "time to reversal (s)") + guides(colour = FALSE, fill = FALSE) + theme(axis.text.x = element_blank()) + scale_y_log10() # glmm was singular glmm <- tyrosine_data %>% lme4::lmer(data = ., log(response.time) ~ food*cond + (1|date)) glmm %>% ref_grid() %>% emmeans::emmeans(pairwise ~ food | cond) lin_mod <- lm(data = tyrosine_data, log(response.time) ~ food*cond) lin_mod %>% ref_grid() %>% emmeans::contrast(method = "pairwise")
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.