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)

#### for 30% octanol ####
SOS_30pct <- readr::read_csv(here::here("extdata","Figure_2D_30pct.csv")) %>%
  dplyr::filter(response.time < 20,
                date != "2019_04_29") %>%
  mutate(plateID = interaction(date, genotype, food, cond, cue))


SOS_30pct %>%
  mutate(food = fct_relevel(food, c("OP50", "JUb39"))) %>%
  dplyr::filter(
  #                     response.time < 20,
  #                     genotype %in% c("tdc-1", "N2"),
                      # food %in% c("OP50", "JUb39")
                      ) %>%
  #                     date %in% c("2019_03_18", "2019_04_09", "2019_04_13", "2019_04_29")) %>%
  ggplot(aes(x = genotype, y = response.time)) +
  #geom_boxplot(outlier.shape =NA) +
  #geom_point(pch = 21) +#, position = position_jitterdodge()) +
  ggbeeswarm::geom_quasirandom(aes(colour = food), width = 0.2) +
  scale_color_plot("grey-blue", drop = TRUE) +
  facet_grid(.~food) +
  add.median('response.time', colour = "red") +
  add.quartiles('response.time') +
  theme_my +
  add.n(genotype) +
  labs(x = "food",
       y = "time to reversal (s)")

library(emmeans)
SOS_30pct %>% 
  dplyr::filter(food != "JUb39; tdcDel::cmR delAADC", date != "2019_04_29") %>% 
  lme4::lmer(data = ., response.time ~ genotype * food + (1|plateID) + (1|date)) %>%
  emmeans(pairwise ~ genotype| food)

library(tidyverse)
library(emmeans)
library(modelr)
library(tidybayes)
library(kableExtra)
library(patchwork)
theme_set(theme_my)

amine_data <- read_csv(here::here('extdata/Figure_2C.csv')) %>%
  mutate(food = fct_relevel(food, "OP50"),
         genotype = fct_relevel(genotype, c("N2", "tdc-1", "tbh-1")))

stan_mod <- rstanarm::stan_glmer(data = amine_data, log(response.time) ~ food * genotype + (food|date) + (1|plateID),
                       family = gaussian,
                       chains = 6,
                       cores = 6)

# fitted <- amine_data %>%
#   data_grid(food, genotype) %>%
#   add_fitted_draws(stan_mod, re_formula = NA) %>%
#   mutate(response.time = exp(.value), data_type = "fit")

fitted1 <- emmeans(stan_mod, pairwise ~ (food | genotype))$contrasts %>% 
  coda::as.mcmc() %>% 
  bayesplot::mcmc_intervals_data(prob = 0.66, prob_outer = 0.95) %>%
  mutate(genotype = levels(amine_data$genotype),
         data_type = "fit",
         food = "JUb39") %>%
  mutate(food = factor(food, levels = c("OP50", "JUb39")),
         genotype = fct_relevel(genotype, c("N2", "tdc-1", "tbh-1"))) %>% 
  mutate_if(is.numeric, funs(. * -1))

sjstats::equi_test(stan_mod)

# frequentist glmm (log transformed)
#singular fit so cannot do full model
lmer <- lme4::lmer(data = amine_data, log(response.time) ~ food * genotype + (1|date) + (1|plateID))

lmer %>% emmeans(pairwise ~ food | genotype)


#plot bars
# plot <- amine_data %>%
#   mutate(response.time = case_when(
#     response.time < 1 ~ 1, 
#     TRUE ~ response.time
#   ),
#   food = fct_relevel(food, "OP50")) %>%
#   ggplot(aes(x = data_type, y = response.time)) +
#   stat_summary(geom = "bar", fun.y = mean, aes(fill = food), width = 0.75, alpha = 0.5) +
#   ggbeeswarm::geom_quasirandom(aes(colour = food), width = 0.2, alpha = 0.75) +
#   scale_color_plot("grey-blue", drop = TRUE) +
#   scale_fill_plot("grey-blue", drop = TRUE) +
#   facet_grid(.~genotype + food, 
#              labeller = labeller(genotype = label_wrap_gen(10), 
#                                  food = as_labeller(c("OP50" = "",
#                                                     "JUb39" = "")))) +
#   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 = "grey") +
#   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 = "grey") +
#   labs(x = "genotype",
#        y = "time to reversal (s)") +
#   guides(colour = FALSE) +
#   theme(axis.text.x = element_blank()) +
#   scale_y_log10()


plot <- amine_data %>%
  format_SOS(., day_correct = genotype) %>%
  ggplot(aes(x = data_type)) +
  geom_relLatency(fitted = fitted1,
                  fillvar = food,
                  dotvar = food,
                  yvar = rel_log) +
  scale_color_plot("grey-blue", drop = TRUE) +
  scale_fill_plot("grey-blue", drop = TRUE) +
  facet_grid(.~genotype+food) +
  add.n('data_type', y.pos = -1.6) +
  labs(x = "genotype",
       y = "relative reversal latency [log(s)]") +
  guides(colour = FALSE) +
  theme(axis.text.x = element_blank())

gt <- ggplot_gtable(ggplot_build(plot))
gt$widths[c(8,12,16,20)] = 5*gt$widths[c(8,12,16,20)]
gt$widths[c(6,10,14,16,18)] = .3*gt$widths[c(6,10,14,16,18)]
grid::grid.draw(gt)

# F-test for interaction term using loci as predictors:
lmer2 <- amine_data %>%
  lme4::lmer(data = ., log(response.time) ~ food * tdc * tbh + (1|date) + (1|plateID))
car::Anova(lmer2, test = "F")


# comaprison of mcmc intervals for supplement:
mcmc.comps <- emmeans(stan_mod, ~ food | genotype, type = "response") %>%
  contrast(method = "pairwise") %>%
  coda::as.mcmc() #%>%
p1 <- bayesplot::mcmc_areas(mcmc.comps)
p2 <- bayesplot::mcmc_intervals(mcmc.comps,prob = 0.66, prob_outer = 0.95)

p1 + p2 + theme(axis.text.y = element_blank())

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