knitr::opts_chunk$set(echo = TRUE) knitr::opts_chunk$set(message = FALSE) knitr::opts_chunk$set(warning = FALSE) library(ProvidenciaChemo) library(rstanarm) library(tidybayes) library(modelr) library(tidyverse) library(emmeans)
library(rstanarm) library(tidybayes) library(modelr) choice <- read_csv(here::here("extdata/Figure_4E_octr1_Choice.csv")) %>% mutate(food = fct_relevel(food, c("OP50","JUb39")), genotype = fct_relevel(genotype, "N2", "octr-1"), N_total = N_OP50 + N_Test, index = (N_Test - N_OP50) / N_total, plateID = factor(seq(1:nrow(.))), p = N_Test / N_total, logit_p = boot::logit(p), strain = food, data_type = "raw")
glmm_mod3 <- lme4::glmer(data = choice %>% filter( N_total > 10, is.na(note)) %>% mutate(strain = fct_relevel(strain, "OP50"), genotype = fct_relevel(genotype, "N2", "octr-1")), cbind(N_Test, N_OP50) ~ strain*genotype + (1|plateID) + (1|date), family = binomial) #%>% summary(glmm_mod3) glmm_mod3 %>% emmeans::emmeans(~ strain | genotype) %>% emmeans::contrast(method = "pairwise") glmm_mod3 %>% emmeans::emmeans(~ genotype | strain) %>% emmeans::contrast(method = "pairwise") #------------bayesian mod for experiment 1------------ stan_mod2 <- rstanarm::stan_glmer(data = choice, cbind(N_Test, N_OP50) ~ strain*genotype + (1|plateID) + (1|date), prior = normal(location = 0.6, scale = 0.6^2), family = binomial, cores = 6, chains = 6, adapt_delta = 0.99) fitted <- choice %>% data_grid(strain, genotype) %>% add_fitted_draws(stan_mod2, re_formula = NA) %>% mutate(logit_p = boot::logit(.value), data_type = "fit") fitted1 <- emmeans(stan_mod2, pairwise ~ (strain | genotype))$contrasts %>% #emmeans(stan_mod2, ~ (strain | genotype)) %>% coda::as.mcmc() %>% bayesplot::mcmc_intervals_data(prob = 0.66, prob_outer = 0.95) %>% mutate(genotype = factor(levels(choice$genotype), levels = levels(choice$genotype)), data_type = "fit", strain = "JUb39") %>% mutate(strain = factor(strain, levels = c("OP50", "JUb39"))) %>% mutate_if(is.numeric, funs(. * -1))
#second plot for tbh-1, octr-1 #generate relative index values: means <- choice %>% filter(strain == "OP50") %>% group_by(genotype, date) %>% summarise(meanOP50 = mean(logit_p)) choice <- full_join(choice, means) %>% mutate(rel.Logit = logit_p - meanOP50) %>% droplevels() plot1 <- choice %>% droplevels() %>% ggplot(aes(x = data_type)) + geom_hline(yintercept = 0, lty = 2, alpha = 0.2) + geom_boxplot(aes(fill = strain, y = rel.Logit), outlier.shape = NA, alpha = 0.75) + ggbeeswarm::geom_quasirandom(aes(colour = strain, y = rel.Logit), width = 0.1, alpha = 0.5) + facet_grid(. ~ genotype + strain) + scale_x_discrete(labels = function(strain) str_wrap(strain, width = 10)) + labs(y = "Relative JUb39 preference (log-odds ratio)") + guides(colour = FALSE) + scale_color_plot(palette = "grey-blue", drop = TRUE) + scale_fill_plot(palette = "grey-blue", drop = TRUE) + # stat_pointinterval(aes(y = logit_p, x = 1.4), # data = fitted, fatten_point = 0, # size_range = c(0.3, 1), colour = "grey" # ) + # stat_summary( # data = fitted, # aes(y = logit_p, x = 1.4), # fun.y = median, # fun.ymin = median, # fun.ymax = median, # geom = "crossbar", # width = 0.05, # alpha = 0.75, # lwd = 0.35, # colour = "darkgrey" # ) + geom_linerange( data = fitted1, aes( x = 1.5, ymin = ll, ymax = hh ), lwd = 0.35, colour = "grey34" ) + geom_crossbar( data = fitted1, aes( y = m, x = 1.5, ymin = l, ymax = h ), width = 0.05, lwd = 0.35, alpha = 0.5, colour = "grey69", fill = "grey69" ) + guides(fill = FALSE) #+ # figure.axes() + #coord_cartesian(ylim = c(-.5, 1)) + #add.n(data_type) space_facets(plot1, n_facets = 3, major_div = 2)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.