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)
glycerol <- read_csv(here::here('extdata/Figure_Ext1b_osm.csv')) %>%
    mutate(food = fct_relevel(food, "OP50"))


#singular fit
lme4::glmer(data = glycerol, cbind(Nout, Ntotal - Nout) ~ food + (1|plate) + (1|date), family = "binomial") %>% summary()

lm(data = glycerol, fraction_out ~ food) %>% summary()

stan_mod <- rstanarm::stan_glmer(data = glycerol, cbind(Nout, Ntotal - Nout) ~ food + (1|plate) + (1|date), family = "binomial",
                    seed = 6459,
                    chains = 4,
                    cores = 4,
                    adapt_delta = 0.99)

fitted <- stan_mod %>% 
  emmeans::emmeans(~food) %>% 
  coda::as.mcmc() %>% 
  bayesplot::mcmc_intervals_data(transformations = function(x) boot::inv.logit(x)) %>%
  mutate(food = factor(c("OP50", "JUb39"), levels = c("OP50", "JUb39")))

glycerol  %>%
  ggplot(aes(y = fraction_out, x = factor(1))) +
  ggbeeswarm::geom_quasirandom(aes(color = food), width = 0.1) +
  scale_color_plot(palette = "2-Ps", drop = TRUE) +
  #theme_black() +
  add.mean(fraction_out, colour = "red", width = 0.1) +
  stat_summary(fun.data = mean_se,
               geom = "errorbar",
               width = 0.2,
               lwd = 0.35,
               colour = "black") +
  facet_grid(~food) +
  guides(colour = FALSE) +
  add.n(1, y.pos = -0.1) +
  coord_cartesian(ylim = c(-.1, 1)) +
  geom_crossbar(data = fitted,
                aes(x = 1.4,
                    y = m,
                    ymin = m,
                    ymax = m),
                width = 0.05,
                colour = "darkgrey") +
  geom_pointrange(data = fitted,
                  aes(x =1.4,
                      y = m,
                      ymin = ll,
                      ymax = hh),
                  size = 0.5,
                  fatten = 0,
                  color = "darkgrey") +
  geom_pointrange(data = fitted,
                  aes(x = 1.4,
                      y = m,
                      ymin = l,
                      ymax =h),
                  size = 1, fatten = 0, 
                  colour = "darkgrey")


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