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(patchwork)
library(tidybayes)
library(modelr)
library(emmeans)
library(magrittr)
library(emmeans)


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

MOb_data <- read.csv(filepath) %>%
  dplyr::filter(date != "5_12_17") %>%
  format_AvoidData(day.correct = "OP50") %>%
  mutate(plate = factor(seq(1:nrow(.)))) %>% droplevels()

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

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

fitted <- recenter_fittedValues(MOb_data, stan_glmm, day.correct = "OP50", BayesFit = "HDI")

MOb_plot <- MOb_data %>%
  plot_plasticityIndex(xvar = strain, 
                       plotColors = "grey-blue", 
                       BayesFit = TRUE,
                       width = 0.2, 
                       alpha = 0.7, 
                       bar = TRUE, 
                       n_pos = -2)  + 
  labs(title = "1-oct (corrected)") + 
  scale_color_plot(palette = "2-Ps", drop = TRUE) +
  scale_fill_plot(palette = "2-Ps", drop = TRUE) +
  guides(color = FALSE) +
  labs(fill="strain") +
  facet_grid(~strain) + 
  coord_cartesian(ylim = c(-2.5, 3.25)) +
  figure.axes()

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


MOb_plot + 
  geom_text(data = glm.contrasts.1, 
            aes(label = paste0("P~",round(p.value, 3)), 
                y = max(1.3*MOb_data$rel.Logit)))

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