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

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

lid_data <- read_csv(filepath) %>%
  format_AvoidData(day.correct = "treatment") %>%
  mutate(plate = factor(seq(1:nrow(.))),
         treatment = factor(treatment, levels = c("none", "lid")),
         data_type = "raw") %>% 
  mutate(group.id = interaction(strain, treatment)) %>% droplevels()

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


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

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

plot1 <- lid_data %>%
  plot_plasticityIndex(xvar = data_type, 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,3))


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



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

plot1

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