knitr::opts_chunk$set(echo = TRUE) knitr::opts_chunk$set(message = FALSE) knitr::opts_chunk$set(warning = FALSE) library(ProvidenciaChemo) library(rstanarm) library(tidybayes) library(modelr) library(tidyverse) library(emmeans)
library(rstanarm) library(tidybayes) library(modelr) TAGOF <- readr::read_csv(here::here("extdata/Figure_4B_octr1TA.csv")) %>% mutate(cond = fct_relevel(cond, "Tyrosine"), TA = case_when( cond == "Tyrosine" ~ 0, cond == "Tyr+10mM_TA" ~ 10 ), data_type = "raw") %>% droplevels() stan_glm <- rstanarm::stan_lmer(TAGOF, formula = log(response.time) ~ genotype*cond + (1 | date) + (1|plateID), #formula = log(response.time) ~ genotype*cond + (1|plateID), seed = 6994, cores = 4, chains = 4, adapt_delta = 0.999) #-----------posterior HDI estimates for relative TA effects--------- TAGOF <- TAGOF %>% # mutate(response.time = case_when( # response.time < 1 ~ 1, # TRUE ~ response.time)) %>% group_by(cond, date) %>% summarise(meanTA_0 = mean(log(response.time))) %>% filter(cond == "Tyrosine") %>% ungroup() %>% select(date, meanTA_0) %>% full_join(., TAGOF) %>% mutate(rel_log = log(response.time) - meanTA_0, food = fct_relevel(food, "OP50")) contrasts.95 <- modelbased::estimate_contrasts(stan_glm, ci = 0.95)[c(2,5),] %>% mutate(genotype = c("N2", "octr-1"), cond = "Tyr+10mM_TA") %>% select(genotype, cond, m = Median, ll = CI_low, hh = CI_high) contrasts.66 <- modelbased::estimate_contrasts(stan_glm, ci = 0.66)[c(2,5),] %>% mutate(genotype = c("N2", "octr-1"), cond = "Tyr+10mM_TA") %>% select(genotype, cond, m = Median, l = CI_low, h = CI_high) fitted <- full_join(contrasts.95, contrasts.66) %>% mutate(genotype = fct_relevel(genotype, "N2"), cond = fct_relevel(cond, "Tyrosine")) %>% mutate_if(is.numeric, function(.) {-1 * .}) #plot with posterior draws fordata normalized to 0mM TA plot <- TAGOF %>% ggplot(aes(x = data_type)) + geom_relLatency(fitted = fitted, fillvar = food, dotvar = food, yvar = rel_log) + scale_fill_manual(values = "#2F8A34") + scale_color_manual(values = "#2F8A34") + facet_grid(.~genotype + cond) + 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 <- TAGOF %>% # data_grid(genotype,cond, date = "2019_05_28", .model = stan_glm) %>% # add_fitted_draws(stan_glm, all_new_levels = TRUE) %>% # mutate(response.time = exp(.value), data_type = "fit") # # plot <- TAGOF %>% # 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, fill = "#2F8A34", width = 0.6, alpha = 0.5) + # ggbeeswarm::geom_quasirandom(colour = "#2F8A34", width = 0.2, alpha = 0.75) + # facet_grid(.~genotype + cond, # labeller = labeller(genotype = label_wrap_gen(10), # food = as_labeller(c("OP50" = "", # "JUb39" = "")))) + # scale_color_plot(palette = "grey-blue-light", drop = TRUE) + # scale_fill_plot(palette = "grey-blue-light", 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() #singular fit TAGOF %>% lme4::lmer(data = ., log(response.time) ~ genotype*cond + (1|plateID) + (1|date)) %>% summary() #singular fit TAGOF %>% lme4::lmer(data = ., log(response.time) ~ genotype*cond + (cond|date)) #%>% #singular fit TAGOF %>% lme4::lmer(data = ., log(response.time) ~ genotype*cond + (1|plateID)) %>% summary() glm_mod <- TAGOF %>% lm(data = ., log(response.time) ~ genotype*cond) glm_mod %>% emmeans::ref_grid() %>% contrast(method = "pairwise")
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.