## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## ----include = FALSE----------------------------------------------------------
# limit threads to avoid CPU time issues on CRAN
data.table::setDTthreads(2)

## ----setup, echo = FALSE------------------------------------------------------
library(designit)
library(tidyverse)

## -----------------------------------------------------------------------------
set.seed(2307111)

conditions <- c("DMSO", sprintf("Compound%02d", 1:11))
# set up batch container
bc <- BatchContainer$new(
  dimensions = list(
    row = 8, col = 12
  )
) |>
  # assign samples with conditions and true effects
  assign_in_order(
    data.frame(
      SampleIndex = 1:96,
      Compound = factor(rep(conditions, 8), levels = conditions),
      trueEffect = rnorm(96, mean = 10, sd = 1)
    )
  )

## -----------------------------------------------------------------------------
# get observations with batch effect
get_observations <- function(bc) {
  bc$get_samples() |>
    mutate(
      plateEffect = 0.5 * sqrt((row - 4.5)^2 + (col - 6.5)^2),
      measurement = trueEffect + plateEffect
    )
}

## -----------------------------------------------------------------------------
dat <- get_observations(bc)

head(dat) |> gt::gt()

## ----rawPlatePlots, fig.height=5.5, fig.width=8-------------------------------
cowplot::plot_grid(
  plotlist = list(
    plot_plate(dat,
      plate = plate,
      row = row, column = col, .color = Compound,
      title = "Layout by treatment"
    ),
    plot_plate(dat,
      plate = plate, row = row, column = col, .color = trueEffect,
      title = "True effect"
    ),
    plot_plate(dat,
      plate = plate, row = row, column = col, .color = plateEffect,
      title = "Plate effect"
    ),
    plot_plate(dat,
      plate = plate, row = row, column = col, .color = measurement,
      title = "Measurement"
    )
  ), ncol = 2, nrow = 2
)

## -----------------------------------------------------------------------------
summary(aov(trueEffect ~ Compound, data = dat))

## -----------------------------------------------------------------------------
summary(aov(measurement ~ Compound, data = dat))

## -----------------------------------------------------------------------------
versusDMSO <- paste0(conditions[-1], "-", conditions[1])
trueDiff <- TukeyHSD(aov(
  trueEffect ~ Compound,
  data = dat
))$Compound
trueDiff[versusDMSO, ]

## -----------------------------------------------------------------------------
measureDiff <- TukeyHSD(aov(measurement ~ Compound,
  data = dat
))$Compound
measureDiff[versusDMSO, ]

## ----boxplot, fig.height=5, fig.width=5---------------------------------------
ggplot(
  dat,
  aes(x = Compound, y = measurement)
) +
  geom_boxplot() +
  ylab("Measurement [w/o randomization]") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1))

## ----eval=FALSE---------------------------------------------------------------
# set.seed(2307111)
# 
# bc_rnd <- optimize_design(
#   bc,
#   scoring = mk_plate_scoring_functions(bc,
#     row = "row", column = "col",
#     group = "Compound"
#   )
# )

## ----include=FALSE------------------------------------------------------------
# this is quite slow, we use cached results
# bc_rnd$get_samples(include_id=TRUE) |> pull(.sample_id) |> dput()

bc_rnd <- bc$move_samples(
  location_assignment =
    c(
      1, 24, 57, 4, 91, 94, 8, 47, 27, 26, 66, 65, 53,
      67, 87, 13, 42, 60, 38, 86, 58, 21, 88, 71, 82, 18,
      56, 11, 77, 64, 31, 45, 85, 25, 3, 36, 69, 75, 50,
      96, 46, 83, 52, 89, 79, 78, 20, 92, 35, 2, 73, 32,
      16, 9, 34, 63, 54, 41, 84, 19, 90, 40, 23, 55, 61,
      29, 12, 68, 74, 39, 70, 33, 80, 5, 48, 15, 93, 49,
      30, 10, 59, 7, 14, 28, 62, 22, 43, 6, 51, 44, 81,
      72, 17, 76, 95, 37
    )
)

## -----------------------------------------------------------------------------
dat_rnd <- get_observations(bc_rnd)

dat_rnd |>
  head() |>
  gt::gt()

## ----randomPlatePlots, fig.height=5.5, fig.width=8----------------------------
cowplot::plot_grid(
  plotlist = list(
    plot_plate(dat_rnd,
      plate = plate,
      row = row, column = col, .color = Compound,
      title = "Layout by treatment"
    ),
    plot_plate(dat_rnd,
      plate = plate, row = row, column = col, .color = trueEffect,
      title = "True effect"
    ),
    plot_plate(dat_rnd,
      plate = plate, row = row, column = col, .color = plateEffect,
      title = "Plate effect"
    ),
    plot_plate(dat_rnd,
      plate = plate, row = row, column = col, .color = measurement,
      title = "Measurement"
    )
  ), ncol = 2, nrow = 2
)

## -----------------------------------------------------------------------------
randMeasureDiff <- TukeyHSD(aov(measurement ~ Compound,
  data = dat_rnd
))$Compound
randMeasureDiff[versusDMSO, ]

## ----randBoxplot, fig.height=5, fig.width=5-----------------------------------
ggplot(
  dat_rnd,
  aes(x = Compound, y = measurement)
) +
  geom_boxplot() +
  ylab("Measurement [with randomization]") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1))