skpr

library(skpr)

#Generate a candidate set of all potential design points to be considered in the experiment
#The hypothetical experiment is determining what affects the caffeine content in coffee
candidate_set = expand.grid(temp = c(80,90,100), 
                            type = c("Kona","Java"),
                            beansize = c("Large","Medium","Small"))

candidate_set
#>    temp type beansize
#> 1    80 Kona    Large
#> 2    90 Kona    Large
#> 3   100 Kona    Large
#> 4    80 Java    Large
#> 5    90 Java    Large
#> 6   100 Java    Large
#> 7    80 Kona   Medium
#> 8    90 Kona   Medium
#> 9   100 Kona   Medium
#> 10   80 Java   Medium
#> 11   90 Java   Medium
#> 12  100 Java   Medium
#> 13   80 Kona    Small
#> 14   90 Kona    Small
#> 15  100 Kona    Small
#> 16   80 Java    Small
#> 17   90 Java    Small
#> 18  100 Java    Small

#Generate the design (default D-optimal)
design = gen_design(candidateset = candidate_set, 
                    model = ~temp + type + beansize,
                    trials=12)

design
#>    temp type beansize
#> 1    80 Kona   Medium
#> 2   100 Java    Small
#> 3    80 Java    Large
#> 4   100 Kona    Large
#> 5   100 Java   Medium
#> 6    80 Kona    Small
#> 7   100 Kona    Small
#> 8    80 Kona    Large
#> 9   100 Java    Large
#> 10   80 Java   Medium
#> 11  100 Kona   Medium
#> 12   80 Java    Small

#Evaluate power for the design with an allowable type-I error of 5%
eval_design(RunMatrix = design,
            model = ~temp + type + beansize,
            alpha=0.05)
#>     parameter            type     power
#> 1 (Intercept)    effect.power 0.8424665
#> 2        temp    effect.power 0.8424665
#> 3        type    effect.power 0.8424665
#> 4    beansize    effect.power 0.5165386
#> 5 (Intercept) parameter.power 0.8424665
#> 6        temp parameter.power 0.8424665
#> 7       type1 parameter.power 0.8424665
#> 8   beansize1 parameter.power 0.5593966
#> 9   beansize2 parameter.power 0.5593966

#Evaluate power for the design using a Monte Carlo simulation. 
#Here, we set the effect size (here, the signal-to-noise ratio) to 1.5.
eval_design_mc(RunMatrix = design,
               model = ~temp + type + beansize,
               alpha=0.05,
               effectsize=1.5)
#>     parameter               type power
#> 1 (Intercept) parameter.power.mc 0.611
#> 2        temp parameter.power.mc 0.623
#> 3       type1 parameter.power.mc 0.625
#> 4   beansize1 parameter.power.mc 0.347
#> 5   beansize2 parameter.power.mc 0.338

#Evaluate power for the design using a Monte Carlo simulation, for a non-normal response. 
#Here, we also increase the number of simululations to improve the precision of the results.
eval_design_mc(RunMatrix = design,
               model = ~temp + type + beansize,
               nsim=5000,
               glmfamily = "poisson",
               alpha=0.05,
               effectsize=c(2,6))
#>     parameter               type  power
#> 1 (Intercept) parameter.power.mc 0.9964
#> 2        temp parameter.power.mc 0.9796
#> 3       type1 parameter.power.mc 0.9766
#> 4   beansize1 parameter.power.mc 0.8854
#> 5   beansize2 parameter.power.mc 0.7088

#skpr was designed to operate with the pipe (%>%) in mind. 
#Here is an example of an entire design of experiments analysis in three lines:

library(dplyr)

expand.grid(temp = c(80,90,100), type = c("Kona","Java"), beansize = c("Large","Medium","Small")) %>%
  gen_design(model = ~temp + type + beansize + beansize:type + I(temp^2), trials=24, optimality="I") %>%
  eval_design_mc(model = ~temp + type + beansize + beansize:type + I(temp^2), alpha=0.05)
#>         parameter               type power
#> 1     (Intercept) parameter.power.mc 0.900
#> 2            temp parameter.power.mc 0.898
#> 3           type1 parameter.power.mc 0.997
#> 4       beansize1 parameter.power.mc 0.917
#> 5       beansize2 parameter.power.mc 0.904
#> 6       I(temp^2) parameter.power.mc 0.636
#> 7 type1:beansize1 parameter.power.mc 0.909
#> 8 type1:beansize2 parameter.power.mc 0.911