Development
In adas.utils: Design of Experiments and Factorial Plans Utilities

library(devtools)
library(usethis)
library(tidyverse)

Imported functions

use_import_from(
  "glue", 
  "glue"
)

use_import_from(
  "readr", c("read_csv", "read_csv2", "write_csv", "write_csv2")
)

use_import_from(
  "tibble", 
  c("tibble", "as_tibble", "rownames_to_column")
)

use_import_from("gghalfnorm", "gghalfnorm")

use_import_from("lubridate", "stamp")

use_import_from(
  "tidyr",
  "pivot_longer"
)

use_import_from(
  "grDevices", c("grey")
)

use_import_from(
  "stats", 
  c("as.formula", "coef", "ecdf", "effects", "na.omit",
               "pnorm", "sd", "setNames", "terms", "aov", "TukeyHSD")
)

use_import_from(
  "utils", c("tail")
)

use_import_from(
  "dplyr", 
   c(
     "across",
     "add_row",
     "arrange", 
     "bind_rows",
     "c_across",
     "desc",
     "filter",
     "if_else",
     "mutate", 
     "n", 
     "pull",
     "relocate",
     "rename",
     "rename_with",
     "rowwise",
     "select",
     "slice_tail",
     "summarize_all",
     "ungroup"
   )
)

use_import_from(
  "purrr", 
   c(
     "accumulate",
     "discard",
     "keep",
     "list_merge",
     "map", 
     "set_names",
     "walk",
     "iwalk"
   )   
)

use_import_from(
  "ggplot2", 
  c(
    "aes", 
    "coord_cartesian",
    "coord_flip",
    "element_text",
    "fortify",
    "ggplot",
    "geom_col",
    "geom_function",
    "geom_hline",
    "geom_label",
    "geom_line",
    "geom_point",
    "geom_qq",
    "geom_qq_line",
    "geom_tile",
    "geom_errorbar",
    "labs",
    "position_dodge",
    "scale_fill_viridis_d",
    "scale_y_continuous",
    "sec_axis",
    "stat_qq",
    "theme"

  )
)

use_import_from(
  "rlang", 
  c("is_formula","sym", ":=")
)

use_import_from(
  "stringr",
  "str_remove_all"
)

Build doc and install

Run this to update documentation and install locally.

devtools::document()
devtools::build_vignettes()
devtools::install()

Examples

library(adas.utils)

fp_design_matrix(5) %>% 
  fp_fraction(~A*B*C*D) %>% 
  fp_fraction(~B*C*D*E) %>% 
  dplyr::mutate(Y=rnorm(dplyr::n()))

fp_alias(~A*B*C*D)
fp_alias(~B*C*D*E)

dm <- fp_design_matrix(3)

fp_augment_center <- function(dm, rep=5) {
  stopifnot("factorial.plan" %in% class(dm))
  r <- nrow(dm)
  fct <- attr(dm, "factors")

  dm %>% 
    add_row(
      StdOrder = (r+1):(r+rep),
      RunOrder = sample((r+1):(r+rep)),
      .treat = "0",
      .rep = 1:rep,
    ) %>% 
    mutate(
      across({fct}, ~ 0)
    )
}

dm %>% 
  fp_augment_center(5)

Making data

$2^2$ CCD

set.seed(0)

f <- function(a, b) {
  1 + 2*a + 3*a^2+ 3*b + 0.05*b^2 + 4*a*b + rnorm(length(a))
}

dm <- fp_design_matrix(2, rep=3) %>% 
  fp_augment_center(rep=4) %>%
  fp_augment_axial(rep=2) %>%
  mutate(
    Y = f(A, B)
  )

dm

dm %>% 
  filter(.treat != "center" & .treat != "axial") %>% 
  lm(Y ~ A*B, data=.) %>%
  anova()

dm %>% 
  filter(.treat != "axial") %>% 
  lm(Y ~ A*I(A^2)*B, data=.) %>%
  anova()

dm %>% 
  lm(Y ~ A*I(A^2)*B*I(B^2)+A:B, data=.) %>%
  anova()

dm %>% 
  lm(Y ~  A * B * I(A^2) * I(B^2), data=.) %>%
  anova()

ccd_experiment_yield <- list(
  base = dm %>% 
    filter(.treat != "center" & .treat != "axial") %>% 
    pull(Y),
  center = dm %>% 
    filter(.treat == "center") %>% 
    pull(Y),
  axial = dm %>%
    filter(.treat == "axial") %>% 
    pull(Y)
)

dm <- fp_design_matrix(3, rep=2)

# fp_add_scale <- function(dm, ..., suffix="_s") {
#   attr(dm, "scales") <- list()
#   for (i in 1:...length()) {
#     name <- ...names()[i]
#     rng <- ...elt(i)
#     if (!(is.numeric(rng) & length(rng) == 2 & is.numeric(dm[[name]]))) {
#       warning("Skipping factor ", name, " (it is not a number, or wrong scale range/type provided)\n")
#       next
#     }
#     dm <- dm %>% 
#       mutate(
#         !!paste0(name, suffix) := scales::rescale(!!sym(name), to=rng)
#       )
#     attr(dm, "scales") <- append(attr(dm, "scales"), setNames(list(rng), name))
#   }
#   return(dm)
# }

dms <-  dm %>% 
  fp_add_scale(A=c(2, 12), B=c(40, 60), suffix="")

dms

fp_design_matrix(2) %>% 
  fp_add_names(A="Temperature", B="Pressure")

dm <- fp_design_matrix(2) %>% 
  fp_add_names(A="Temperature", B="Pressure") %>% 
  fp_add_scale(A=c(2, 12), B=c(40, 60), suffix="_s") %>%
  fp_write_csv("design_matrix.csv")

dm %>%
  fp_read_csv("design_matrix.csv")

Centered design

set.seed(0)
fp <- fp_design_matrix(2, rep=3) %>% 
  mutate(Y=f(A, B))
fp

fp %>% 
  lm(Y ~ A*B, data=.) %>% 
  anova()

All factors and their interactions are significant. But is the two-level model enough? Let's check for the quadratic terms, by augmenting the plan with a central point repeated 4 times. We also load the center field from the ccd_experiment_yield dataset:

set.seed(0)
f <- function(a, b) {
  1 + 2*a + 3*b + (3*a^2 + 0.05*b^2)*0.5 + 4*a*b + rnorm(length(a))
}

fpc <- fp %>% 
  fp_augment_center(rep=4) %>% 
  mutate(Y=f(A,B))

fp <- fpc %>% 
  filter(.treat != "center")

fpc

fpc %>% 
  lm(Y ~ A*B+I(A^2), data=.) %>% 
  anova()

fp %>% 
  lm(Y~A*B, data=.) %>% 
  predict(newdata=fpc, interval="confidence") %>% 
  bind_cols(fpc) %>% 
  filter(.treat == "center") %>% 
  summarise(lwr=min(lwr), upr=max(upr)) %>% 
  mutate(what="base") %>% 
  bind_rows(
    fpc %>% 
      filter(.treat == "center") %>% 
      pull(Y) %>% 
      t.test() %>% 
      broom::tidy() %>% 
      select(lwr=conf.low, upr=conf.high) %>% 
      mutate("what"="center")
  ) %>% 
  ggplot(aes(x=what, ymin=lwr, ymax=upr)) +
  geom_errorbar()