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)))
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.