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(emmeans)
library(modelr)
library(tidybayes)
library(kableExtra)
library(patchwork)
theme_set(theme_my)

receptors_exp1 <- read_csv(here::here('extdata/Figure_Ext9b_receptors.csv')) %>%
  filter(response.time < 20,
         genotype %in% c("N2", "tyra-2", "ser-3", "octr-1")) %>%
  mutate(food = fct_relevel(food, "OP50"),
         genotype = fct_relevel(genotype, c("N2", "octr-1", "ser-3", "tyra-2")),
         data_type = "raw") %>% droplevels()  

stan_mod <- rstanarm::stan_glmer(data = filter(receptors_exp1) %>% droplevels(),
                                 log(response.time) ~ food * genotype  + (1|plateID), #add food|date back
                       family = gaussian,
                       seed = 695,
                       chains = 6,
                       cores = 6)

## ---------- Bayesian credible intervals 
fitted <- receptors_exp1 %>%
  data_grid(food, genotype) %>%
  add_fitted_draws(stan_mod, re_formula = NA) %>%
  mutate(response.time = exp(.value), data_type = "fit")

fitted1 <- emmeans(stan_mod, pairwise ~ (food | genotype))$contrasts %>% 
  coda::as.mcmc() %>% 
  bayesplot::mcmc_intervals_data(prob = 0.66, prob_outer = 0.95) %>%
  mutate(genotype = factor(levels(receptors_exp1$genotype), levels = levels(receptors_exp1$genotype)),
         data_type = "fit",
         food = "JUb39") %>%
  mutate(food = factor(food, levels = c("OP50", "JUb39"))) %>% 
  mutate_if(is.numeric, funs(. * -1))

#--------relative effect plot (log)

plot <- receptors_exp1 %>%
  format_SOS(., day_correct = genotype) %>%
  ggplot(aes(x = data_type)) +
  #geom_hline(yintercept = 0, linetype = "dashed", colour = "grey") +
  geom_relLatency(fitted = fitted1,
                  fillvar = food,
                  dotvar = food,
                  yvar = rel_log) +
  scale_color_plot("grey-blue", drop = TRUE) +
  scale_fill_plot("grey-blue", drop = TRUE) +
  facet_grid(.~genotype+food) +
  add.n('data_type', y.pos = -1.6) +
  labs(x = "genotype",
       y = "relative reversal latency [log(s)]") +
  guides(colour = FALSE) +
  theme(axis.text.x = element_blank())

space_facets(plot, n_facets = 4)

sjstats::equivalence_test(stan_mod)

# frequentist glmm (log transformed) - singular fit
lmer <- lme4::lmer(data = receptors_exp1, log(response.time) ~ food * genotype + (food|date) + (1|plateID))

# frequentist glmm (log transformed) non-singular

lmer<- lme4::lmer(data = receptors_exp1, log(response.time) ~ food * genotype + (1|date) + (1|plateID))
# lmer_octr <- lme4::lmer(data = receptors_exp1 %>% mutate(genotype = fct_relevel(genotype, "octr-1")), log(response.time) ~ food * genotype + (1|date) + (1|plateID))

#lmer_octr %>% emmeans(pairwise ~ food | genotype)
lmer %>% emmeans(pairwise ~ genotype | food)
lmer %>% emmeans(pairwise ~ food | genotype)

# t-test for predictors using Kenward-Rogers estimated degrees of freedom
lmer %>% lmerTest::as_lmerModLmerTest() %>% summary(ddf = "Kenward-Roger")

# parametric bootstrap to estimate octr-1 interaction effect:
library(lme4)
lmer_add <- lme4::lmer(data = receptors_exp1, log(response.time) ~ food + genotype + (1|date) + (1|plateID))
lmer_boot <- lme4::bootMer(x=lmer, FUN=fixef, re.form = NA, nsim=200)
boot::boot.ci(lmer_boot, index=6, type = "perc", conf = 0.95)

# comparison of mcmc intervals for supplement:
mcmc.comps <- emmeans(stan_mod, ~ food | genotype, type = "response") %>%
  contrast(method = "pairwise") %>%
  coda::as.mcmc() #%>%
p1 <- bayesplot::mcmc_areas(mcmc.comps)
p2 <- bayesplot::mcmc_intervals(mcmc.comps,prob = 0.66, prob_outer = 0.95)

p1 + p2 + theme(axis.text.y = element_blank())

lmer %>% emmeans::ref_grid() %>% emmeans::contrast(method = "pairwise", adjust = "fdr") %>% kableExtra::kable() %>% kable_styling(bootstrap_options = c("striped", "hover"))

#test interaction effects:
#phia::testInteractions(lmer_octr, fixed = "genotype", across = "food")

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