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(emmeans)
library(modelr)
library(tidybayes)
library(kableExtra)
library(patchwork)
theme_set(theme_my)

library(tidyverse)
library(patchwork)
library(emmeans)
library(tidybayes)

filepath <- here::here("extdata/Figure_1C_antibiotic.csv")

Gent_data <- read_csv(filepath) %>%
  mutate(
    treatment = case_when(
           treatment == "live" ~ "none",
           treatment == "dead" ~ "dead"
         ),
    col_map = interaction(treatment, strain)) %>% droplevels() %>%
  format_AvoidData(day.correct = "treatment") %>%
  mutate(plate = factor(seq(1:nrow(.))),
        treatment = fct_relevel(treatment, "none"),
         food = fct_relevel(strain, "OP50"),
         data_type = "raw") %>%
  filter(!is.na(date))


glmm <- lme4::glmer(data = Gent_data, 
                    formula = cbind(nCue,nControl) ~ treatment * food + (1|date) + (1|plate),
                    family = binomial,
                    control= lme4::glmerControl(optimizer="bobyqa"))

stan_mod <-  rstanarm::stan_glmer(data = Gent_data, 
                  formula = cbind(nCue,nControl) ~ food * treatment + (1 | date) + (1|plate),
                  chains = 6, cores = 4, seed = 2000,iter=6000,
                  family = binomial,
    control = list(adapt_delta=0.99))

fitted <- recenter_fittedValues(Gent_data, stan_mod, BayesFit = "HDI", day.correct = "treatment")

plot1 <- Gent_data %>%
    ggplot(aes(x = data_type, y = rel.Logit)) +
    stat_summary(geom = "bar", fun.y = mean, aes(fill = food), width = 0.5, alpha = 0.75) +
    # add.mean(logit_p, colour = "red") +
    # add.quartiles(logit_p) +
    ggbeeswarm::geom_quasirandom(aes(colour = food), width = 0.1, alpha = 0.75) +
    stat_summary(geom = "errorbar", fun.data = mean_se, width = 0.4) +
    facet_grid(.~ treatment + food, scales = "free_x") +
    scale_x_discrete(labels = function(food) str_wrap(food, width = 10)) +
    labs(y = "Test bacteria preference (log-odds)") +
    #theme(panel.spacing = unit(4, "lines")) +
  #stat_summary(fun.data = mean_se, geom = "errorbar", width = 0.2) +
  guides(colour = FALSE) +
  scale_color_plot(palette = "grey-blue-green", drop = TRUE) +
    scale_fill_plot(palette = "grey-blue-green", drop = TRUE) +
    stat_pointinterval(aes(y=rel.Logit, x = 1.4),
                       data = fitted, fatten_point = 0,
                       size_range = c(0.3, 1), colour = "grey") +
    stat_summary(data = fitted,
                 aes(y=rel.Logit, x = 1.4),
                 fun.y = median,
                 fun.ymin = median,
                 fun.ymax = median,
                 geom = "crossbar",
                 width = 0.05,
                 lwd = 0.35,
                 colour = "grey") +
    guides(fill = FALSE)


plot1 = Gent_data %>%
  plot_plasticityIndex(xvar = treatment, BayesFit = TRUE, width = 0.2, alpha = 0.7, bar = TRUE, n_pos = -2, dot_color = col_map) +
  facet_grid(.~food + treatment) +
  labs(title = "1-oct (corrected)") +
  scale_color_plot(palette = "2-each", drop = TRUE, rev = TRUE) +
  scale_fill_plot(palette = "2-each", drop = TRUE, rev = TRUE) +
  guides(color = FALSE, fill = FALSE) + coord_cartesian(ylim = c(-2.5, 3.25))


glm.lsm.list <- emmeans::ref_grid(glmm) %>% emmeans::emmeans(., pairwise ~ treatment | food, adjust = "mvt")
glm.contrasts <- update(glm.lsm.list$contrast, adjust = "mvt", by.vars = NULL)

library(nlme)
lm <- lm(rel.Logit ~ food * treatment,data = Gent_data)
anova(lm)
lm.lsm.list <- emmeans::ref_grid(lm) %>% emmeans::emmeans(., pairwise ~ treatment | food, adjust = "Tukey")

plot1

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