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)
library(rstanarm)
theme_set(theme_my)

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

supe_data <- read_csv(filepath) %>%
  filter(assay == "oct", feed_state == "fed") %>%
  format_AvoidData(day.correct = "treatment", min_p = 0.01) %>%
  mutate(plate = factor(seq(1:nrow(.))),
         #treatment = factor(treatment, levels = c("none", "lid")),
         data_type = "raw") %>% 
  mutate(group.id = interaction(strain, treatment),
         food = fct_relevel(strain, "OP50"),
         treatment = fct_relevel(treatment, "none", "NGM", "OPCM"),
         feed_state = fct_relevel(feed_state, "fed")) %>% droplevels()

#convergence errors with date random effect, left out of this model
glmm <- lme4::glmer(data = supe_data, 
                    formula = cbind(nCue,nControl) ~ strain * treatment + (1|plate) + (1|date),
                    family = binomial,
                    control= lme4::glmerControl(optimizer="bobyqa"))


stan_glmm <-  rstanarm::stan_glmer(data = supe_data,
                  formula = cbind(nCue,nControl) ~ strain * treatment + (1|plate) + (1|date),
                  chains = 4, cores = 4, seed = 2000,iter=6000,
                  family = binomial,
                  prior_intercept = normal(0.005, .002),
    control = list(adapt_delta=0.99))

fitted <- recenter_fittedValues(supe_data, stan_glmm, day.correct = "treatment", BayesFit = "fitted_draws")

plot1 <- supe_data %>%
  filter(assay == "oct", feed_state == "fed") %>%
  plot_plasticityIndex(xvar = strain, width = 0.2, BayesFit = TRUE, bar = TRUE, alpha = 0.75) +
  scale_color_plot(palette = "grey-blue", drop = TRUE) + 
  scale_fill_plot(palette = "grey-blue", drop = TRUE) +
  guides(color = FALSE) + 
  facet_grid(~treatment + strain) + coord_cartesian(ylim = c(-2,4))


#glm.contrasts <- emmeans::emmeans(glmm, pairwise ~ strain | treatment)
emmeans::ref_grid(glmm) %>% emmeans::contrast(., method = "pairwise") %>%
  broom::tidy() #%>% kable() %>% kable_styling()


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

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