In mikeod38/ProvidenciaChemo:

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")

mikeod38/ProvidenciaChemo documentation built on April 6, 2020, 11:57 p.m.