Tooth Growth

rm(list = ls())
library(ggpubr)
library(rstatix)
library(grid)
library(gridExtra)
require(bayesrules)
require(rstanarm)
require(bayesplot)
require(tidybayes)
library(brms)
library(modelr)
library(MCMCglmm)
library(ggstatsplot)
theme_set(theme_tidybayes() + panel_border())

df <- ToothGrowth
df$dose <- as.factor(df$dose)
head(df, 10)
p <- ggboxplot(ToothGrowth, 
               x = "supp", 
               y = "len", 
               color = "supp", 
               palette = "jco", 
               add = "jitter")
dev.off()
ggbetweenstats(
  data  = df,
  x     = supp,
  y     = len,
  title = "Distribution of len across supp")

grouped_ggbetweenstats(
  data             = df,
  x                = dose,
  y                = len,
  grouping.var     = supp,
  outlier.tagging  = TRUE,
  outlier.coef     = 2,
  ggsignif.args    = list(textsize = 4, tip_length = 0.01),
  p.adjust.method  = "bonferroni",
  palette          = "default_jama",
  package          = "ggsci",
  plotgrid.args    = list(nrow = 1),
  annotation.args  = list(title = "Differences in len by dose for different supp"))

ggwithinstats(
  data    = df,
  x       = dose,
  y       = len,
  title   = "Tooth growth")

grouped_ggwithinstats(
  data            = df,
  x               = dose,
  y               = len,
  type            = "np",
  xlab            = "Dose",
  ylab            = "Length",
  grouping.var    = supp,
  outlier.tagging = TRUE)

gghistostats(
  data       = df,
  x          = len,
  title      = "Length",
  test.value = 15,
  binwidth   = 1)

grouped_gghistostats(
  data              = df,
  x                 = len,
  test.value        = 15,
  type              = "nonparametric",
  xlab              = "Length",
  grouping.var      = supp,
  normal.curve      = TRUE,
  normal.curve.args = list(color = "red", size = 1),
  ggtheme           = ggthemes::theme_tufte(),
  plotgrid.args     = list(nrow = 1),
  annotation.args   = list(title = "Length by supp"))


p1 <- ggwithinstats(
  data             = df,
  x                = dose,
  y                = len,
  type = "p",
  effsize.type = "d",
  conf.level = 0.99,
  title = "Parametric test",
  package = "ggsci",
  palette = "nrc_npg")
p2 <- ggwithinstats(
  data             = df,
  x                = dose,
  y                = len,
  xlab            = "Dose",
  ylab            = "Length",
  type = "np",
  conf.level = 0.99,
  title = "Non-parametric Test",
  package = "ggsci",
  palette = "uniform_startrek")
p3 <- ggwithinstats(
  data             = df,
  x                = dose,
  y                = len,
  xlab            = "Dose",
  ylab            = "Length",
  type = "r",
  conf.level = 0.99,
  title = "Robust Test",
  package = "wesanderson",
  palette = "Royal2")
p4 <- ggwithinstats(
  data             = df,
  x                = dose,
  y                = len,
  xlab            = "Dose",
  ylab            = "Length",
  type = "bayes",
  title = "Bayesian Test",
  package = "ggsci",
  palette = "nrc_npg")
combine_plots(
  plotlist = list(p1, p2, p3, p4),
  plotgrid.args = list(nrow = 2),
  annotation.args = list(
    title = "Effect of disgust on dose on length ",
    caption = "Source: Toothgrow dataset"))



p <- ggboxplot(ToothGrowth, x = "supp", y = "len",
               color = "supp", palette = "jco",
               add = "jitter")
#  Add p-value
p + stat_compare_means()
# Change method
g1<-p + stat_compare_means(method = "t.test") + theme_bw()
g2<-p + stat_compare_means(method = "anova") + theme_bw()
g3<-p + stat_compare_means(method = "wilcox.test") + theme_bw()
grid.arrange(g1,g2,g3,ncol=3)

compare_means(len ~ supp, 
              data = ToothGrowth, 
              paired = FALSE)

g<-ggplot(ToothGrowth, aes(
  x = len, 
  fill = as.factor(supp)))+
  geom_density(alpha=0.5)+
  labs(fill="supp")+
  theme_bw()
set_palette(g,palette = "jco")
ggpaired(ToothGrowth, x = "supp", y = "len",
         color = "supp", line.color = "gray", line.size = 0.4,
         palette = "jco")+
  stat_compare_means(paired = TRUE) + theme_bw()
compare_means(len ~ supp, 
              data = ToothGrowth, 
              paired = TRUE)

my_comparisons <- list( c("0.5", "1"), c("1", "2"), c("0.5", "2") )
ggboxplot(ToothGrowth, x = "dose", y = "len",
          color = "dose", palette = "jco")+ 
  stat_compare_means(comparisons = my_comparisons, label.y = c(29, 35, 40))+
  stat_compare_means(label.y = 45) +
  theme_bw()
g<-ggplot(ToothGrowth, aes(
  x = len, 
  fill = as.factor(dose)))+
  geom_density(alpha=0.5)+
  labs(fill="dose")+
  theme_bw()
set_palette(g,palette = "jco")

stat.test <- df %>%
  group_by(dose) %>%
  t_test(len ~ supp) %>%
  adjust_pvalue() %>%
  add_significance("p.adj")
stat.test <- stat.test %>% 
  add_xy_position(x = "supp")
bxp <- ggboxplot(df, x = "supp", y = "len", fill = "supp",
                 facet.by = "dose")
p<-bxp + 
  stat_pvalue_manual(stat.test, label = "p.adj") +
  scale_y_continuous(expand = expansion(mult = c(0.05, 0.10)))
set_palette(p,palette = "jco")

g<-ggplot(ToothGrowth, aes(
  x = len, 
  fill = as.factor(supp)))+
  geom_density(alpha=0.5)+
  facet_grid(~dose)+
  labs(fill="supp")+
  theme_bw()
set_palette(g,palette = "jco")

# Statistical test
stat.test <- df %>%
  t_test(len ~ supp, paired = TRUE) %>%
  add_significance()
bxp <-  ggpaired(df, x = "supp", y = "len", fill = "#E7B800",
                 line.color = "gray", line.size = 0.4)
stat.test <- stat.test %>% add_xy_position(x = "supp")
bxp + stat_pvalue_manual(stat.test, label = "p.signif")
bxp + 
  stat_pvalue_manual(stat.test, label = "{p}{p.signif}") +
  scale_y_continuous(expand = expansion(mult = c(0.05, 0.10)))

## De Bayesiaanse manier
df <- ToothGrowth
df$dose<-as.numeric(df$dose);str(df)
get_prior(len ~ dose + supp, 
          data = df, 
          family = gaussian)
df_normal_prior<-rstanarm::stan_glm(len ~ dose + supp, 
                    data = df, 
                    family = gaussian,
                    prior_intercept = normal(7, 1.5),
                    prior = normal(0, 2.5, autoscale = TRUE),
                    prior_aux = rstanarm::exponential(1, autoscale = TRUE),
                    chains = 4, iter = 5000*2, seed = 84735, 
                    prior_PD = TRUE)
prior_summary(df_normal_prior)
plot(df_normal_prior)
df_normal <- update(df_normal_prior,
                    prior_PD = FALSE)
plot(df_normal, plotfun="dens")
df_normal
pp_check(df_normal, 
         plotfun = "hist", 
         nreps = 15) + 
  geom_vline(xintercept = 0)
pp_check(df_normal, 
         plotfun = "dens", 
         nreps = 15) + 
  geom_vline(xintercept = 0)

mcmc_trace(df_normal)
mcmc_dens_overlay(df_normal)
mcmc_acf(df_normal)

df %>%
  add_fitted_draws(df_normal, n = 50) %>%
  ggplot(aes(x = dose, y = len, color = supp)) +
  geom_line(aes(y = .value, group = paste(supp, .draw)), 
            alpha = .2) +
  geom_point(data = df, size = 0.5) + theme_bw()

mn_prediction <- posterior_predict(df_normal, 
                                   newdata = data.frame(dose = 1.8,
                                                        supp = "OJ"))
mn_prediction 
mcmc_hist(mn_prediction, binwidth = 1) + 
  geom_vline(xintercept = 0, lty=2, color="red") + 
  xlab("Predicted number of len for dose = 1.8 and supp = OJ")

mn_prediction <- posterior_predict(df_normal, 
                                   newdata = data.frame(dose = c(1.8,1.8),
                                                        supp = c("OJ","VC")), 
                                   draws=1000)
mn_prediction<-as.data.frame(mn_prediction)
mn_prediction$diff<-mn_prediction$`1`-mn_prediction$`2`
melt_mn_pred<-reshape2::melt(mn_prediction);str(melt_mn_pred)
ggplot(melt_mn_pred, aes(x=value, 
                         fill=as.factor(variable)))+
  geom_density(alpha=0.5)+
  labs(fill="supp", x="len", title="Predicted len for dose = 1.8")+
  theme_bw()
set.seed(84735)
normal_predictions <- posterior_predict(df_normal, newdata = df)
ppc_intervals_grouped(df$len, yrep = normal_predictions, 
                      x = df$dose, 
                      group = df$supp,
                      prob = 0.5, prob_outer = 0.95,
                      facet_args = list(scales = "fixed"))
prediction_summary(model = df_normal, 
                   data = df)
set.seed(84735)
normal_cv <- prediction_summary_cv(model = df_normal, 
                                   data = df, k = 10)
normal_cv$cv


df$dose_fc<-as.factor(df$dose);str(df)
priorget <- get_prior(len ~ dose_fc + supp, 
                      data = df, family = gaussian)  
priorget
df_normal_prior2 = brm(
  len ~ dose_fc + supp, 
  data = df, 
  prior = c(
    prior(normal(7, 1.5), class=Intercept),
    prior(normal(0, 2.5), class=b, coef=dose_fc1),
    prior(normal(0, 2.5), class=b, coef=dose_fc2),
    prior(normal(0, 2.5), class=b, coef=suppVC),
    prior(exponential(1), class=sigma)),
  control = list(adapt_delta = .99))
priorget
df_normal_prior2$prior
summary(df_normal_prior2)
plot(df_normal_prior2)
ndraws = 4000
p = df %>%
  group_by(supp) %>%
  add_epred_draws(df_normal_prior2, ndraws = ndraws)%>%
  ggplot(., aes(x=.epred, fill=supp))+
  geom_density(alpha=0.5)+
  facet_wrap(~dose_fc, scales = "free")+
  labs(x="Predicted len", 
       title="Distributions of len by supp and dose provided via Bayesian inference (ndraws=4000)")+
  theme_bw()
set_palette(p,palette = "jco")

get_prior(len ~ dose + supp, 
          data = df, family=gaussian())
df_normal_prior3 = brm(
  len ~ dose + supp, 
  data = df,
  prior = c(
    prior(normal(7, 1.5), class=Intercept),
    prior(normal(4, 1), class=b, coef=dose),
    prior(normal(3, 0.5), class=b, coef=suppVC),
    prior(exponential(1), class=sigma)),
  control = list(adapt_delta = .99))
get_prior(len ~ dose + supp, 
          data = df, family=gaussian())
df_normal_prior3$prior
get_variables(df_normal_prior3)
summary(df_normal_prior3)
ndraws = 4000
try<-df_normal_prior3 %>%
  recover_types(df)%>%
  spread_draws(b_dose, b_suppVC)
ggplot(try, aes(x=b_suppVC, fill=as.factor(.chain)))+
  geom_density(alpha=0.5)+
  labs(fill="Chain")+
theme_bw()
ggplot(try, aes(x=b_dose, 
                y=b_suppVC, 
                color=as.factor(.chain)))+
  geom_point(alpha=0.7)+
  labs(color="Chain")+
  theme_bw()

df <- ToothGrowth
df$dose<-as.numeric(df$dose);str(df)
priorget <- get_prior(len ~ (1|supp) + (1|dose), 
                      data = df, family = gaussian)  
priorget
df$len
df_normal_prior4 = brm(
  len ~ (1|supp) + (1|dose), 
  data = df, 
  prior = c(
    prior(normal(7, 1.5), class=Intercept),
    prior(gamma(3,3), class=sd, coef=Intercept, group=dose),
    prior(exponential(2), class=sd, coef=Intercept, group=supp),
    prior(exponential(1), class=sigma)),
  control = list(adapt_delta = .99))
priorget
df_normal_prior4$prior
summary(df_normal_prior4)

df %>%
  recover_types(df) %>%
  add_epred_draws(df_normal_prior4, 
                  dpar = c("mu", "sigma")) %>%
  ggplot(aes(x = supp)) +
  stat_dist_slab(aes(dist = "norm", arg1 = mu, arg2 = sigma), 
                 slab_color = "gray65", alpha = 1/10, fill = NA, 
                 data = . %>% sample_draws(100), scale = .5) +
  stat_halfeye(aes(y = .epred), side = "bottom", scale = .5) +
  geom_point(aes(y = len), data = df, shape = 21, fill = "#9ECAE1", size = 2, position = position_nudge(x = -.2))+
  facet_wrap(~dose)+
  theme_bw()

df_normal_prior4 %>%
  spread_draws(r_supp[supp]) %>%
  compare_levels(r_supp, by = supp) %>%
  dplyr::ungroup() %>%
  mutate(supp = reorder(supp, r_supp)) %>%
  ggplot(aes(y = supp, x = r_supp)) +
  stat_halfeye() +
  geom_vline(xintercept = 0, linetype = "dashed") +theme_bw()

df_normal_prior4 %>%
  spread_draws(r_dose[dose]) %>%
  compare_levels(r_dose, by = dose, comparison="ordered") %>%
  dplyr::ungroup() %>%
  mutate(dose = reorder(dose, r_dose)) %>%
  ggplot(aes(y = dose, x = r_dose)) +
  stat_halfeye() +
  geom_vline(xintercept = 0, linetype = "dashed") +theme_bw()

df_normal_prior4 %>%
  spread_draws(r_dose[dose]) %>%
  compare_levels(r_dose, by = dose, comparison="pairwise") %>%
  dplyr::ungroup() %>%
  mutate(dose = reorder(dose, r_dose)) %>%
  ggplot(aes(y = dose, x = r_dose)) +
  stat_halfeye() +
  geom_vline(xintercept = 0, linetype = "dashed") +theme_bw()

df_normal_prior4 %>%
  spread_draws(r_dose[dose]) %>%
  compare_levels(r_dose, by = dose, comparison="control") %>%
  dplyr::ungroup() %>%
  mutate(dose = reorder(dose, r_dose)) %>%
  ggplot(aes(y = dose, x = r_dose)) +
  stat_halfeye() +
  geom_vline(xintercept = 0, linetype = "dashed") +theme_bw()

get_variables(df_normal_prior4)
summary(df_normal_prior4)
ToothGrowth %>%
  add_predicted_draws(df_normal_prior4) %>%
  ggplot(aes(x = len, y = .prediction, 
             color = ordered(dose), 
             fill = ordered(dose))) +
  stat_lineribbon(
                  .width = c(.95, .80, .50), 
                  alpha = 1/4) +
  geom_point(color="black", alpha=0.2) +
  scale_fill_brewer(palette = "Set2") +
  scale_color_brewer(palette = "Dark2")+
  theme_bw()+
  facet_wrap(~dose, scales="free")

ToothGrowth %>%
  add_predicted_draws(df_normal_prior4) %>%
  ggplot(aes(x = len, y = .prediction, 
             color = ordered(supp), 
             fill = ordered(supp))) +
  stat_lineribbon(
    .width = c(.95, .80, .50), 
    alpha = 1/4) +
  geom_point(color="black", alpha=0.2) +
  scale_fill_brewer(palette = "Set2") +
  scale_color_brewer(palette = "Dark2")+
  theme_bw()+
  facet_wrap(~supp, scales="free")