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)

choice <- read_csv(here::here("extdata/Figure_4E_octr1_Choice.csv")) %>%
  mutate(food = fct_relevel(food, c("OP50","JUb39")),
         genotype = fct_relevel(genotype, "N2", "octr-1"),
         N_total = N_OP50 + N_Test,
    index = (N_Test - N_OP50) / N_total, 
    plateID = factor(seq(1:nrow(.))),
    p = N_Test / N_total,
    logit_p = boot::logit(p),
    strain = food,
    data_type = "raw")

glmm_mod3 <- lme4::glmer(data = choice %>%
                           filter(
                                  N_total > 10,
                                  is.na(note)) %>%
                           mutate(strain = fct_relevel(strain, "OP50"),
                                  genotype = fct_relevel(genotype, "N2", "octr-1")),
                         cbind(N_Test, N_OP50) ~ strain*genotype + (1|plateID) + (1|date), family = binomial) #%>%

summary(glmm_mod3)

glmm_mod3 %>% emmeans::emmeans(~ strain | genotype) %>% emmeans::contrast(method = "pairwise")

glmm_mod3 %>% emmeans::emmeans(~ genotype | strain) %>% emmeans::contrast(method = "pairwise")

#------------bayesian mod for experiment 1------------
stan_mod2 <- rstanarm::stan_glmer(data = choice,
                     cbind(N_Test, N_OP50) ~ strain*genotype + (1|plateID) + (1|date),
                     prior = normal(location = 0.6, scale = 0.6^2),
                     family = binomial,
                     cores = 6,
                     chains = 6,
                     adapt_delta = 0.99)

fitted <- choice %>% 
  data_grid(strain, genotype) %>%
  add_fitted_draws(stan_mod2, re_formula = NA) %>%
  mutate(logit_p = boot::logit(.value), data_type = "fit")


fitted1 <- emmeans(stan_mod2, pairwise ~ (strain | genotype))$contrasts %>% 
  #emmeans(stan_mod2,  ~ (strain | genotype)) %>%
  coda::as.mcmc() %>% 
  bayesplot::mcmc_intervals_data(prob = 0.66, prob_outer = 0.95) %>%
  mutate(genotype = factor(levels(choice$genotype), levels = levels(choice$genotype)),
         data_type = "fit",
         strain = "JUb39") %>%
  mutate(strain = factor(strain, levels = c("OP50", "JUb39"))) %>% 
  mutate_if(is.numeric, funs(. * -1))

#second plot for tbh-1, octr-1

#generate relative index values:
means <- choice %>%
      filter(strain == "OP50") %>%
      group_by(genotype, date) %>%
      summarise(meanOP50 = mean(logit_p))

choice <- full_join(choice, means) %>% mutate(rel.Logit = logit_p - meanOP50) %>% droplevels()


plot1 <- choice %>%
  droplevels() %>%
  ggplot(aes(x = data_type)) +
  geom_hline(yintercept = 0, lty = 2, alpha = 0.2) +
  geom_boxplot(aes(fill = strain, y = rel.Logit), outlier.shape = NA, alpha = 0.75) +
  ggbeeswarm::geom_quasirandom(aes(colour =  strain, y = rel.Logit), width = 0.1, alpha = 0.5) +
  facet_grid(. ~ genotype + strain) +
  scale_x_discrete(labels = function(strain) str_wrap(strain, width = 10)) +
  labs(y = "Relative JUb39 preference (log-odds ratio)") +
  guides(colour = FALSE) +
  scale_color_plot(palette = "grey-blue", drop = TRUE) +
  scale_fill_plot(palette = "grey-blue", drop = TRUE) +
  # stat_pointinterval(aes(y = logit_p, x = 1.4),
  #   data = fitted, fatten_point = 0,
  #   size_range = c(0.3, 1), colour = "grey"
  # ) +
  # stat_summary(
  #   data = fitted,
  #   aes(y = logit_p, x = 1.4),
  #   fun.y = median,
  #   fun.ymin = median,
  #   fun.ymax = median,
  #   geom = "crossbar",
  #   width = 0.05,
  #   alpha = 0.75,
  #   lwd = 0.35,
  #   colour = "darkgrey"
  #  ) +
  geom_linerange(
    data = fitted1,
    aes(
      x = 1.5,
      ymin = ll,
      ymax = hh
    ),
    lwd = 0.35,
    colour = "grey34"
  ) +
  geom_crossbar(
    data = fitted1,
    aes(
      y = m,
      x = 1.5,
      ymin = l,
      ymax = h
    ),
    width = 0.05,
    lwd = 0.35,
    alpha = 0.5,
    colour = "grey69",
    fill = "grey69"
  ) +
  guides(fill = FALSE) #+
  # figure.axes() +
  #coord_cartesian(ylim = c(-.5, 1)) +
  #add.n(data_type)

space_facets(plot1, n_facets = 3, major_div = 2)

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