In mikeod38/ProvidenciaChemo:

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

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