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