## ----setup, include = FALSE, warning = FALSE----------------------------------
knitr::opts_chunk$set(comment = FALSE, 
                      warning = FALSE, 
                      message = FALSE)

## -----------------------------------------------------------------------------
library(cellmig)
library(ggplot2)
library(ggforce)
ggplot2::theme_set(new = theme_bw(base_size = 10))

## ----eval=FALSE---------------------------------------------------------------
# if (!require("BiocManager", quietly = TRUE))
#     install.packages("BiocManager")
# 
# BiocManager::install("cellmig")

## -----------------------------------------------------------------------------
data("d", package = "cellmig")
str(d)

## ----fig.width=7, fig.height=6------------------------------------------------
ggplot(data = d)+
  facet_wrap(facets = ~paste0("compound=", compound), 
             scales = "free_y", ncol = 2)+
  geom_sina(aes(x = as.factor(dose), col = plate, y = v, group = well), 
            size = 0.5)+
  theme_bw()+
  theme(legend.position = "top",
        strip.text.x = element_text(margin = margin(0.03,0,0.03,0, "cm")))+
  ylab(label = "migration velocity")+
  xlab(label = '')+
  scale_color_grey()+
  guides(color = guide_legend(override.aes = list(size = 3)))+
  guides(shape = guide_legend(override.aes = list(size = 3)))+
  scale_y_log10()+
  annotation_logticks(base = 10, sides = "l")

## ----fig.width=7, fig.height=6------------------------------------------------
dm <- aggregate(v~well+plate+compound+dose, data = d, FUN = mean)
ggplot(data = dm)+
  facet_wrap(facets = ~paste0("compound=", compound), 
             scales = "free_y", ncol = 2)+
  geom_sina(aes(x = as.factor(dose), col = plate, y = v, group = well), 
            size = 1.5, alpha = 0.7)+
  theme_bw()+
  theme(legend.position = "top",
        strip.text.x = element_text(margin = margin(0.03,0,0.03,0, "cm")))+
  ylab(label = "migration velocity")+
  xlab(label = '')+
  scale_color_grey()+
  guides(color = guide_legend(override.aes = list(size = 3)))+
  guides(shape = guide_legend(override.aes = list(size = 3)))+
  scale_y_log10()+
  annotation_logticks(base = 10, sides = "l")

## ----fig.width=7, fig.height=3.5----------------------------------------------
o <- cellmig(x = d,
             control = list(mcmc_warmup = 300, # nr. of MCMC warmup step?
                            mcmc_steps = 1000, # nr. of MCMC iteration steps?
                            mcmc_chains = 2,   # nr. of MCMC chains
                            mcmc_cores = 2))   # nr. of MCMC cores

## -----------------------------------------------------------------------------
str(o$posteriors$delta_t)

## ----fig.width=6, fig.height=3.3----------------------------------------------
ggplot(data = o$posteriors$delta_t)+
  geom_line(aes(x = dose, y = mean, col = compound, group = compound))+
  geom_point(aes(x = dose, y = mean, col = compound))+
  geom_errorbar(aes(x = dose, y = mean, ymin = X2.5., ymax = X97.5., 
                    col = compound), width = 0.1)+
  ylab(label = expression("Overall treatment effect ("*delta*")"))+
  theme(legend.position = "top")

## ----fig.width=9.5, fig.height=5----------------------------------------------
get_dose_response_profile(x = o)+
  patchwork::plot_layout(widths = c(.7, 1, 4))

## ----fig.width=14, fig.height=6-----------------------------------------------
u <- get_pairs(x = o, exponentiate = FALSE)
u$plot

## ----fig.width=7, fig.height=2.5----------------------------------------------
get_groups(x = o)
u <- get_violins(x = o, 
                 from_groups = get_groups(x = o)$group,
                 to_group = "C2|D1",
                 exponentiate = FALSE)
u$plot

## ----fig.width=5, fig.height=6------------------------------------------------
g <- get_ppc_violins(x = o, wrap = TRUE, ncol = 3)
g+scale_y_log10()

## ----fig.width=4, fig.height=4------------------------------------------------
g <- get_ppc_means(x = o)
g

## ----fig.height=2, fig.width=6------------------------------------------------
g_alpha_p <- ggplot(data = o$posteriors$alpha_p)+
  geom_errorbarh(aes(y = plate_id, x = mean, xmin = X2.5., xmax = X97.5.),
                 height = 0.1)+
  geom_point(aes(y = plate_id, x = mean))

g_sigma <- ggplot()+
  geom_errorbarh(data = o$posteriors$sigma_bio,
                 aes(y = "sigma_bio",
                     x = mean, xmin = X2.5., xmax = X97.5.),
                 height = 0.1)+
  geom_errorbarh(data = o$posteriors$sigma_tech,
                 aes(y = "sigma_tech",
                     x = mean, xmin = X2.5., xmax = X97.5.),
                 height = 0.1)+
  geom_point(data = o$posteriors$sigma_bio,
             aes(y = "sigma_bio", x = mean))+
  geom_point(data = o$posteriors$sigma_tech,
             aes(y = "sigma_tech", x = mean))+
  ylab(label = '')

g_alpha_p|g_sigma

## -----------------------------------------------------------------------------
sessionInfo()