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

## ----echo=FALSE---------------------------------------------------------------
template_dir <- system.file("input_templates/cell_assay_data",
                            package = "RMeDPower2")

## ----eval=FALSE---------------------------------------------------------------
# install.packages("devtools")
# library(devtools)
# install_github('gladstone-institutes/RMeDPower2', build_vignettes=TRUE, force = TRUE)

## ----results='hide'-----------------------------------------------------------
suppressPackageStartupMessages(library(RMeDPower2))
suppressPackageStartupMessages(library(ggplot2))
suppressPackageStartupMessages(library(dplyr))
suppressPackageStartupMessages(library(magrittr))

data <- plate_assay_pilot_data

## -----------------------------------------------------------------------------
knitr::kable(data[1:5,] )

## -----------------------------------------------------------------------------
design <- readDesign(file.path(template_dir,"design_cell_assay.json"))

## -----------------------------------------------------------------------------
model <- readProbabilityModel(file.path(template_dir,"prob_model.json"))

## -----------------------------------------------------------------------------
power_param <- readPowerParams(file.path(template_dir,"power_param.json"))

## ----eval=FALSE---------------------------------------------------------------
# diagnose_res <- diagnoseDataModel(data, design, model)
# data_update = diagnose_res$Data_updated

## ----eval=FALSE---------------------------------------------------------------
# diagnose_res_update <- diagnoseDataModel(data_update, design_update, model_update)

## ----eval=FALSE---------------------------------------------------------------
# power_res <- calculatePower(data, design, model, power_param)

## ----eval=FALSE---------------------------------------------------------------
# estimate_res <- getEstimatesOfInterest(data, design, model, print_plots = F)

## ----echo=FALSE---------------------------------------------------------------
suppressPackageStartupMessages(require(dplyr))
suppressPackageStartupMessages(library(RMeDPower2))
template_dir <- system.file("input_templates/cell_assay_data",
                            package = "RMeDPower2")

data <- plate_assay_pilot_data %>%
  select(experiment, line, classification, cell_size2) %>%
  mutate(experiment = as.factor(experiment), line = as.factor(line), classification = as.factor(classification))

str(data)

## -----------------------------------------------------------------------------
design <- new("RMeDesign")
print(design)

## ----echo=FALSE---------------------------------------------------------------
template_dir <- system.file("input_templates/cell_assay_data",
                            package = "RMeDPower2")

## -----------------------------------------------------------------------------
design <- readDesign(file.path(template_dir, "design_cell_assay.json"))
print(design)

## -----------------------------------------------------------------------------
model <- readProbabilityModel(file.path(template_dir, "prob_model.json"))
print(model)

## -----------------------------------------------------------------------------
power_param <- readPowerParams(file.path(template_dir, "power_param.json"))
print(power_param)

## ----flowchart, echo=FALSE, fig.align = "center", fig.cap="Flowchart illustrating the functionality of the package", fig.small=TRUE, echo=FALSE----
knitr::include_graphics(system.file("figures/package_structure.png", package = "RMeDPower2"))

## ----applications, echo=FALSE, fig.align = "center", fig.cap="Flowchart illustrating the functionality of the package", fig.small=TRUE, echo=FALSE----
knitr::include_graphics(system.file("figures/application_domains.png", package = "RMeDPower2"))

## -----------------------------------------------------------------------------
template_dir <- system.file("input_templates/cell_assay_data", package = "RMeDPower2")
data <- plate_assay_pilot_data
design <- readDesign(file.path(template_dir,"design_cell_assay.json"))
model <- readProbabilityModel(file.path(template_dir,"prob_model.json"))
power_param <- readPowerParams(file.path(template_dir,"power_param.json"))

## ----echo=FALSE---------------------------------------------------------------
print(ggplot(data, aes(x=experiment, y=cell_size2, color = line)) +
  geom_boxplot() +
  theme_bw() +
  theme(text = element_text(size = 16)))

print(ggplot(data, aes(x=cell_size2, color = experiment)) +
  stat_ecdf() +
  theme_bw() +
  theme(text = element_text(size = 16)))

## -----------------------------------------------------------------------------
data %<>% filter(cell_size2 < quantile(cell_size2, 0.95))

## -----------------------------------------------------------------------------
diagnose_res <- diagnoseDataModel(data, design, model)

## -----------------------------------------------------------------------------
print(diagnose_res$models)

## ----output-plots1, results='asis'--------------------------------------------
plot_no <- 1
plots <- diagnose_res$diagnostic_plots$natural_scale$plots
captions <- diagnose_res$diagnostic_plots$natural_scale$captions
for (i in seq_along(plots)) {
  print(plots[[i]])
  cat(sprintf("\n\n\n**Figure %d:** %s. %s\n\n", plot_no,names(plots)[i], captions[i]))
  plot_no <- plot_no + 1
}

cat("\n\n\n")
cooks_plots <- diagnose_res$cooks_plots$natural_scale
for(i in 1:length(cooks_plots)) {
  print(cooks_plots[[i]])
  cat(sprintf("\n\n\n**Figure %d:**\n\n", plot_no))
  plot_no <- plot_no + 1
}

print("Inferred outlier groups are")
print(diagnose_res$inferred_outlier_groups$natural_scale)

## ----output-plots2, results='asis'--------------------------------------------
plots <- diagnose_res$diagnostic_plots$log_scale$plots
captions <- diagnose_res$diagnostic_plots$log_scale$captions
for (i in seq_along(plots)) {
  print(plots[[i]])
  cat(sprintf("\n\n\n**Figure %d:** %s. %s\n\n", plot_no,names(plots)[i], captions[i]))
  plot_no <- plot_no + 1
}


cooks_plots <- diagnose_res$cooks_plots$log_scale

cat("\n\n\n")
for(i in 1:length(cooks_plots)) {
  print(cooks_plots[[i]])
  cat(sprintf("\n\n\n**Figure %d:**\n\n", plot_no))
  plot_no <- plot_no + 1
}

print("Inferred outlier groups are")
print(diagnose_res$inferred_outlier_groups$natural_scale)

## -----------------------------------------------------------------------------
design_update <- design
data_update <- diagnose_res$Data_updated
print(colnames(data_update))
design_update@response_column <-  "cell_size2_logTransformed"
estimate_res <- getEstimatesOfInterest(data_update, design_update, model,print_plots = FALSE)

## ----results="asis", echo = FALSE---------------------------------------------
print(estimate_res[[2]]$plots)
knitr::asis_output(sprintf("\n\n\n**Figure %d:**\n\n", plot_no))
plot_no <- plot_no + 1


## -----------------------------------------------------------------------------
##print model summary output
print(estimate_res[[1]])

## -----------------------------------------------------------------------------
power_param@max_size <- 15
power_res <- calculatePower(data_update, design_update, model, power_param)

## ----results="asis", echo = FALSE---------------------------------------------
print(power_res$plots)
knitr::asis_output(sprintf("\n\n\n**Figure %d:**\n\n", plot_no))
plot_no <- plot_no + 1

## -----------------------------------------------------------------------------
full_data <- plate_assay_full_data

##filter the top 5% of the observations
full_data %<>% filter(perim_2th_effect < quantile(perim_2th_effect, 0.95))

design@response_column <- "perim_2th_effect"
design@condition_column <- "classif"
full_diag_res <- diagnoseDataModel(full_data, design, model)
full_data_update <- full_diag_res$Data_updated
full_design <- design
full_design@response_column <- "perim_2th_effect_logTransformed"

## ----results='asis'-----------------------------------------------------------
cooks_plots <- full_diag_res$cooks_plots$log_scale
cat("\n\n\n")
for(i in 1:length(cooks_plots)) {
  print(cooks_plots[[i]])
  cat(sprintf("\n\n\n**Figure %d:**\n\n", plot_no))
  plot_no <- plot_no + 1
}

print("Inferred outlier groups are")
print(diagnose_res$inferred_outlier_groups$natural_scale)


## -----------------------------------------------------------------------------
full_res <- getEstimatesOfInterest(full_data_update, full_design, model, print_plots = FALSE)
print(full_res[[1]])

## ----results="asis", echo = FALSE---------------------------------------------
print(full_res[[2]]$plots)
knitr::asis_output(sprintf("\n\n\n**Figure %d:**\n\n", plot_no))
plot_no <- plot_no + 1
cat("\n\n\n")

## ----sessionInfo, echo=FALSE--------------------------------------------------
sessionInfo()