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)
theme_set(theme_classic())

library(tidyverse)
library(emmeans)
library(modelr)
library(tidybayes)

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

amine_data <- read_csv(filepath) %>%
  dplyr::filter(strain %in% c("OP50", "JUb39"),
                genotype %in% c("N2","tdc-1","tbh-1", "cat-2", "tph-1")) %>% 
  droplevels() %>%
  format_AvoidData(day.correct = "genotype") %>%
  mutate(plate = factor(seq(1:nrow(.))),
         genotype = factor(genotype, levels = c("N2","tdc-1","tbh-1", "cat-2", "tph-1"))) %>% droplevels()

#### 
glmm <- lme4::glmer(data = amine_data, 
                    formula = cbind(nCue,nControl) ~ genotype * strain + (1 | date) +  (1|plate),
                    family = binomial)

####
glm<- lme4::lmer(data = amine_data, 
                    formula = CI ~ genotype * strain + (1|date))

stan_glmm <-  rstanarm::stan_glmer(data = amine_data,
                  formula = cbind(nCue,nControl) ~
                    genotype * strain + (1 + strain + genotype | date) +  (1|plate),
                  chains = 6, cores = 6, seed = 2000,iter=6000,
                  family = binomial,
    control = list(adapt_delta=0.99))

fitted <- recenter_fittedValues(amine_data, stan_glmm, BayesFit = "fitted_draws", day.correct = "OP50_by_genotype")

# glm.contrasts <- emmeans::ref_grid(glmm) %>% emmeans::contrast(., method = "pairwise") %>%
#   broom::tidy() %>%
#   filter(level1 == c("OP50,none ")) %>%
#   mutate(strain = forcats::as_factor(c("JUb39", "OP50", "JUb39")),
#          treatment = forcats::as_factor(c("none", "lid", "lid")),
#          data_type = forcats::as_factor(rep("raw", 3)))

plot <- plot_plasticityIndex(df = amine_data, xvar = strain, BayesFit = TRUE, width = 0.2, alpha = 0.7, bar = TRUE) + 
  labs(title = "day corrected") +
  facet_grid(.~genotype+strain) + 
  scale_color_plot(palette = "grey-blue", drop = TRUE) + 
  guides(color = FALSE) +
  scale_fill_plot(palette = "grey-blue", drop = TRUE) + 
  guides(fill = FALSE) +
  coord_cartesian(ylim = c(-1.5,3.5)) +
  geom_hline(yintercept = 0, alpha = 0.5)

# plot2 <- amine_data %>%
#   ggplot(aes(x = strain, y = rel.Logit)) +
#   geom_boxplot(aes(fill = strain), alpha = 0.5, width = 0.8, outlier.shape = NA) +
#   ggbeeswarm::geom_quasirandom(aes(colour = strain), width = 0.2, alpha = 0.75) +
#   facet_grid(.~genotype) +
#   scale_fill_plot(palette = "grey-blue", drop = TRUE) +
#   scale_color_plot(palette = "grey-blue", drop = TRUE) +
#   theme(panel.spacing.x = unit(2, "lines"))

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

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

lsm.list <- emmeans::ref_grid(glmm) %>% emmeans::lsmeans(., pairwise ~ strain | genotype) #%>% summary(adjust = "bon")
contrasts <- update(lsm.list$contrast, adjust = "mvt", by.vars = NULL)

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