Fully-crossed simulation designs

Author

Phil Chalmers

In situations where researchers are unaware which variables or combinations of variables influence a given outcome, and when computing the relevant models is very expensive, Monte Carlo simulations can be used. In this case computer-driven experimental designs can be constructed with the goal of detecting (or screening) which factors affect a given outcome variable through statistical inference. Montgomery (2017) describes this class of screening experiments as “… tests in which many factors are considered and the objective is to identify those factors that have significant effects”, which are typically relevant in early stages of a given project.

The most common designs implemented (and indeed, the default used by ) are those which fully-cross the factor variables, leading to a set of \(C_1\times C_2 \times \cdots \times C_n\) unique conditions to study. Fully-crossed designs are beneficial when detecting main effects, two-way interactions, and up to \(n\)-way interactions (provided that the number of replications per unique condition is greater than 1).

For example, consider the following fully-crossed experimental design that combines five simulation factors in the context of an independent samples \(t\)-test’s ability to detect true mean differences: sample size, the magnitude of the mean difference, group variances, distributions, whether the group sizes were equal or unequal, if the distributions were symmetric/asymmetric.

Fully-crossed design example

library(SimDesign)

Design <- createDesign(sample_size=c(100,200),
                       mean_diff=c(.25, 2),
                       variance.ratio=c(1,4),
                       equal_size=c(TRUE, FALSE),
                       dists=c('norm', 'skew'),
                       same_dists=c(TRUE, FALSE),
                       symmetric=c(TRUE, FALSE),
                       # remove same-normal combo
                       subset = !(symmetric & dists == 'norm'))
Design
# A tibble: 96 × 7
   sample_size mean_diff variance.ratio equal_size dists same_dists symmetric
         <dbl>     <dbl>          <dbl> <lgl>      <chr> <lgl>      <lgl>    
 1         100      0.25              1 TRUE       skew  TRUE       TRUE     
 2         200      0.25              1 TRUE       skew  TRUE       TRUE     
 3         100      2                 1 TRUE       skew  TRUE       TRUE     
 4         200      2                 1 TRUE       skew  TRUE       TRUE     
 5         100      0.25              4 TRUE       skew  TRUE       TRUE     
 6         200      0.25              4 TRUE       skew  TRUE       TRUE     
 7         100      2                 4 TRUE       skew  TRUE       TRUE     
 8         200      2                 4 TRUE       skew  TRUE       TRUE     
 9         100      0.25              1 FALSE      skew  TRUE       TRUE     
10         200      0.25              1 FALSE      skew  TRUE       TRUE     
# ℹ 86 more rows
dim(Design)
[1] 96  7
#-------------------------------------------------------------------

Generate <- function(condition, fixed_objects = NULL) {
    Attach(condition)
    N1 <- ifelse(equal_size, sample_size/2, sample_size/4)
    N2 <- sample_size - N1
    DV1 <- if(dists == 'norm') rnorm(N1, sd=1) else rchisq(N1, df=5) - 5
    if(dists == 'norm' && same_dists){ 
        DV2 <- rnorm(N2, sd=sqrt(variance.ratio))
    } else {
        DV2 <- (rchisq(N2, df=5) - 5) * sqrt(variance.ratio)
        if(!symmetric) DV2 <- abs(0 - DV2)
    }
    dat <- data.frame(DV=c(DV1, DV2 + mean_diff), 
                      group=c(rep('g1', N1), rep('g2', N2)))
    dat
}

Analyse <- function(condition, dat, fixed_objects = NULL) {
    test <- t.test(DV ~ group, data=dat, var.equal=TRUE)
    ret <- c(detect=test$p.value < .05)
    ret
}

#-------------------------------------------------------------------

# note the omission of the summarise function
res <- runSimulation(design=Design, replications=20, 
                     generate=Generate, analyse=Analyse, verbose=FALSE)
res
# A tibble: 1,920 × 8
   sample_size mean_diff variance.ratio equal_size dists same_dists symmetric
         <dbl>     <dbl>          <dbl> <lgl>      <chr> <lgl>      <lgl>    
 1         100      0.25              1 TRUE       skew  TRUE       TRUE     
 2         100      0.25              1 TRUE       skew  TRUE       TRUE     
 3         100      0.25              1 TRUE       skew  TRUE       TRUE     
 4         100      0.25              1 TRUE       skew  TRUE       TRUE     
 5         100      0.25              1 TRUE       skew  TRUE       TRUE     
 6         100      0.25              1 TRUE       skew  TRUE       TRUE     
 7         100      0.25              1 TRUE       skew  TRUE       TRUE     
 8         100      0.25              1 TRUE       skew  TRUE       TRUE     
 9         100      0.25              1 TRUE       skew  TRUE       TRUE     
10         100      0.25              1 TRUE       skew  TRUE       TRUE     
# ℹ 1,910 more rows
# ℹ 1 more variable: detect <lgl>

Given these collected results, ANOVA-based approaches can be employed to investigate main and interaction effects for the recorded outcome variable.

# main effects only
mod <- aov(detect ~ ., data=res)
summary(mod)
                 Df Sum Sq Mean Sq F value  Pr(>F)    
sample_size       1    0.2     0.2    2.42   0.120    
mean_diff         1   62.7    62.7  661.28 < 2e-16 ***
variance.ratio    1    3.3     3.3   34.28 5.6e-09 ***
equal_size        1    0.6     0.6    5.98   0.015 *  
dists             1    4.6     4.6   48.57 4.4e-12 ***
same_dists        1    6.2     6.2   65.25 1.2e-15 ***
symmetric         1  107.5   107.5 1133.88 < 2e-16 ***
Residuals      1912  181.3     0.1                    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# all two-way effects
mod2 <- aov(detect ~ (.)^2, data=res)
summary(mod2)
                             Df Sum Sq Mean Sq F value  Pr(>F)    
sample_size                   1    0.2     0.2    3.87 0.04929 *  
mean_diff                     1   62.7    62.7 1056.73 < 2e-16 ***
variance.ratio                1    3.3     3.3   54.77 2.0e-13 ***
equal_size                    1    0.6     0.6    9.56 0.00202 ** 
dists                         1    4.6     4.6   77.62 < 2e-16 ***
same_dists                    1    6.2     6.2  104.27 < 2e-16 ***
symmetric                     1  107.5   107.5 1811.94 < 2e-16 ***
sample_size:mean_diff         1    0.0     0.0    0.01 0.92537    
sample_size:variance.ratio    1    0.0     0.0    0.22 0.63955    
sample_size:equal_size        1    0.0     0.0    0.08 0.77871    
sample_size:dists             1    0.0     0.0    0.36 0.55112    
sample_size:same_dists        1    0.1     0.1    1.48 0.22343    
sample_size:symmetric         1    0.1     0.1    1.07 0.30191    
mean_diff:variance.ratio      1    0.0     0.0    0.43 0.51205    
mean_diff:equal_size          1    0.6     0.6    9.56 0.00202 ** 
mean_diff:dists               1    0.0     0.0    0.74 0.38926    
mean_diff:same_dists          1    7.9     7.9  132.77 < 2e-16 ***
mean_diff:symmetric           1   39.6    39.6  666.44 < 2e-16 ***
variance.ratio:equal_size     1    0.3     0.3    4.64 0.03131 *  
variance.ratio:dists          1    0.1     0.1    1.58 0.20833    
variance.ratio:same_dists     1    0.1     0.1    1.97 0.16012    
variance.ratio:symmetric      1    2.9     2.9   48.98 3.6e-12 ***
equal_size:dists              1    0.2     0.2    3.20 0.07385 .  
equal_size:same_dists         1    0.0     0.0    0.08 0.77871    
equal_size:symmetric          1    0.7     0.7   11.07 0.00089 ***
dists:same_dists              1   16.4    16.4  276.45 < 2e-16 ***
same_dists:symmetric          1    0.1     0.1    2.22 0.13598    
Residuals                  1892  112.3     0.1                    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# all three-way effects
mod3 <- aov(detect ~ (.)^3, data=res)
summary(mod3)
                                        Df Sum Sq Mean Sq F value  Pr(>F)    
sample_size                              1    0.2     0.2    4.66 0.03094 *  
mean_diff                                1   62.7    62.7 1273.28 < 2e-16 ***
variance.ratio                           1    3.3     3.3   66.00 8.1e-16 ***
equal_size                               1    0.6     0.6   11.52 0.00070 ***
dists                                    1    4.6     4.6   93.53 < 2e-16 ***
same_dists                               1    6.2     6.2  125.64 < 2e-16 ***
symmetric                                1  107.5   107.5 2183.25 < 2e-16 ***
sample_size:mean_diff                    1    0.0     0.0    0.01 0.91811    
sample_size:variance.ratio               1    0.0     0.0    0.26 0.60720    
sample_size:equal_size                   1    0.0     0.0    0.10 0.75774    
sample_size:dists                        1    0.0     0.0    0.43 0.51292    
sample_size:same_dists                   1    0.1     0.1    1.79 0.18144    
sample_size:symmetric                    1    0.1     0.1    1.28 0.25715    
mean_diff:variance.ratio                 1    0.0     0.0    0.52 0.47172    
mean_diff:equal_size                     1    0.6     0.6   11.52 0.00070 ***
mean_diff:dists                          1    0.0     0.0    0.89 0.34464    
mean_diff:same_dists                     1    7.9     7.9  159.98 < 2e-16 ***
mean_diff:symmetric                      1   39.6    39.6  803.01 < 2e-16 ***
variance.ratio:equal_size                1    0.3     0.3    5.59 0.01812 *  
variance.ratio:dists                     1    0.1     0.1    1.91 0.16727    
variance.ratio:same_dists                1    0.1     0.1    2.38 0.12312    
variance.ratio:symmetric                 1    2.9     2.9   59.02 2.5e-14 ***
equal_size:dists                         1    0.2     0.2    3.85 0.04976 *  
equal_size:same_dists                    1    0.0     0.0    0.10 0.75774    
equal_size:symmetric                     1    0.7     0.7   13.34 0.00027 ***
dists:same_dists                         1   16.4    16.4  333.11 < 2e-16 ***
same_dists:symmetric                     1    0.1     0.1    2.68 0.10174    
sample_size:mean_diff:variance.ratio     1    0.1     0.1    2.38 0.12312    
sample_size:mean_diff:equal_size         1    0.0     0.0    0.01 0.91811    
sample_size:mean_diff:dists              1    0.3     0.3    5.08 0.02430 *  
sample_size:mean_diff:same_dists         1    0.1     0.1    2.38 0.12312    
sample_size:mean_diff:symmetric          1    0.1     0.1    2.68 0.10174    
sample_size:variance.ratio:equal_size    1    0.0     0.0    0.86 0.35483    
sample_size:variance.ratio:dists         1    0.0     0.0    0.89 0.34464    
sample_size:variance.ratio:same_dists    1    0.0     0.0    0.26 0.60720    
sample_size:variance.ratio:symmetric     1    0.0     0.0    0.14 0.70560    
sample_size:equal_size:dists             1    0.1     0.1    2.33 0.12693    
sample_size:equal_size:same_dists        1    0.0     0.0    0.86 0.35483    
sample_size:equal_size:symmetric         1    0.0     0.0    0.14 0.70560    
sample_size:dists:same_dists             1    0.1     0.1    1.53 0.21656    
sample_size:same_dists:symmetric         1    0.0     0.0    0.02 0.89979    
mean_diff:variance.ratio:equal_size      1    0.0     0.0    0.52 0.47172    
mean_diff:variance.ratio:dists           1    1.2     1.2   23.73 1.2e-06 ***
mean_diff:variance.ratio:same_dists      1    0.9     0.9   17.78 2.6e-05 ***
mean_diff:variance.ratio:symmetric       1    0.5     0.5    9.91 0.00167 ** 
mean_diff:equal_size:dists               1    0.4     0.4    8.04 0.00462 ** 
mean_diff:equal_size:same_dists          1    0.1     0.1    1.28 0.25813    
mean_diff:equal_size:symmetric           1    1.1     1.1   21.71 3.4e-06 ***
mean_diff:dists:same_dists               1   14.6    14.6  296.98 < 2e-16 ***
mean_diff:same_dists:symmetric           1    0.0     0.0    0.02 0.89979    
variance.ratio:equal_size:dists          1    0.0     0.0    0.13 0.71622    
variance.ratio:equal_size:same_dists     1    0.0     0.0    0.26 0.60720    
variance.ratio:equal_size:symmetric      1    0.3     0.3    5.73 0.01681 *  
variance.ratio:dists:same_dists          1    0.5     0.5   10.71 0.00109 ** 
variance.ratio:same_dists:symmetric      1    0.0     0.0    0.14 0.70560    
equal_size:dists:same_dists              1    0.0     0.0    0.43 0.51292    
equal_size:same_dists:symmetric          1    0.0     0.0    0.78 0.37810    
Residuals                             1862   91.7     0.0                    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

and so on.

Fractional factorial experiments

When interest is primarily in detecting main and two-way interaction effects only then substantial gains can be made by using fractional factorial designs. A fractional factorial design is a reduced version of the fully factorial design, allowing for a more efficient use of resources as it reduces the sample size and number of conditions studied by way of discarding higher-order factor combinations.

Note that to use the FrF2 package each factor must consist of exactly two-levels.

library(FrF2)
# help(FrF2)

# 7 factors with resolution 5
fr <- FrF2(nfactors=7, resolution=5)
fr[1:6,]
   A  B  C  D  E  F  G
1  1 -1 -1  1  1  1  1
2  1 -1 -1  1 -1 -1  1
3  1 -1  1  1  1 -1  1
4  1  1  1  1  1  1  1
5 -1  1 -1 -1 -1 -1 -1
6  1  1  1 -1 -1 -1 -1
# Create working simulation design given -1/1 combinations
fDesign <- createDesign(sample_size=c(100,200),
                        mean_diff=c(.25, 2),
                        variance.ratio=c(1,4),
                        equal_size=c(TRUE, FALSE),
                        dists=c('norm', 'skew'),
                        same_dists=c(TRUE, FALSE),
                        symmetric=c(TRUE, FALSE),
                        # remove same-normal combo
                        subset = !(symmetric & dists == 'norm'),
                        fractional=fr)
fDesign
# A tibble: 48 × 7
   sample_size mean_diff variance.ratio equal_size dists same_dists symmetric
         <dbl>     <dbl>          <dbl> <lgl>      <chr> <lgl>      <lgl>    
 1         100      0.25              1 TRUE       skew  FALSE      TRUE     
 2         200      2                 1 FALSE      skew  FALSE      FALSE    
 3         100      2                 4 FALSE      skew  FALSE      FALSE    
 4         100      2                 1 TRUE       skew  TRUE       TRUE     
 5         100      0.25              4 FALSE      norm  FALSE      FALSE    
 6         200      0.25              1 TRUE       norm  TRUE       FALSE    
 7         200      2                 4 FALSE      skew  TRUE       FALSE    
 8         200      2                 4 TRUE       skew  TRUE       TRUE     
 9         200      2                 4 FALSE      norm  FALSE      FALSE    
10         100      0.25              4 FALSE      skew  TRUE       FALSE    
# ℹ 38 more rows
dim(fDesign)
[1] 48  7

Notice that the fractional factor design has fewer condition than the fully-crossed version above. However, despite this reduced number of conditions valid main and two-way interaction effects are still estimable.

# note the omission of the summarise function
fres <- runSimulation(design=fDesign, replications=20, 
                     generate=Generate, analyse=Analyse, verbose=FALSE)
fres
# A tibble: 960 × 8
   sample_size mean_diff variance.ratio equal_size dists same_dists symmetric
         <dbl>     <dbl>          <dbl> <lgl>      <chr> <lgl>      <lgl>    
 1         100      0.25              1 TRUE       skew  FALSE      TRUE     
 2         100      0.25              1 TRUE       skew  FALSE      TRUE     
 3         100      0.25              1 TRUE       skew  FALSE      TRUE     
 4         100      0.25              1 TRUE       skew  FALSE      TRUE     
 5         100      0.25              1 TRUE       skew  FALSE      TRUE     
 6         100      0.25              1 TRUE       skew  FALSE      TRUE     
 7         100      0.25              1 TRUE       skew  FALSE      TRUE     
 8         100      0.25              1 TRUE       skew  FALSE      TRUE     
 9         100      0.25              1 TRUE       skew  FALSE      TRUE     
10         100      0.25              1 TRUE       skew  FALSE      TRUE     
# ℹ 950 more rows
# ℹ 1 more variable: detect <lgl>
# main effects only
mod <- aov(detect ~ ., data=fres)
summary(mod)
                Df Sum Sq Mean Sq F value  Pr(>F)    
sample_size      1    1.4     1.4   14.88 0.00012 ***
mean_diff        1   35.3    35.3  388.81 < 2e-16 ***
variance.ratio   1    1.1     1.1   11.76 0.00063 ***
equal_size       1    0.3     0.3    3.72 0.05403 .  
dists            1    1.9     1.9   20.67 6.2e-06 ***
same_dists       1    4.8     4.8   53.10 6.6e-13 ***
symmetric        1   57.6    57.6  635.03 < 2e-16 ***
Residuals      952   86.4     0.1                    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# all two-way effects
mod2 <- aov(detect ~ (.)^2, data=fres)
summary(mod2)
                            Df Sum Sq Mean Sq F value  Pr(>F)    
sample_size                  1    1.4     1.4   27.23 2.2e-07 ***
mean_diff                    1   35.3    35.3  711.31 < 2e-16 ***
variance.ratio               1    1.1     1.1   21.51 4.0e-06 ***
equal_size                   1    0.3     0.3    6.81  0.0092 ** 
dists                        1    1.9     1.9   37.82 1.2e-09 ***
same_dists                   1    4.8     4.8   97.15 < 2e-16 ***
symmetric                    1   57.6    57.6 1161.76 < 2e-16 ***
sample_size:mean_diff        1    0.3     0.3    6.81  0.0092 ** 
sample_size:variance.ratio   1    0.1     0.1    2.10  0.1475    
sample_size:equal_size       1    0.1     0.1    1.34  0.2465    
sample_size:dists            1    0.3     0.3    6.05  0.0141 *  
sample_size:same_dists       1    0.0     0.0    0.76  0.3847    
sample_size:symmetric        1    1.6     1.6   32.27 1.8e-08 ***
mean_diff:variance.ratio     1    0.2     0.2    4.12  0.0427 *  
mean_diff:equal_size         1    0.0     0.0    0.34  0.5622    
mean_diff:dists              1    0.1     0.1    2.06  0.1516    
mean_diff:same_dists         1    5.1     5.1  102.95 < 2e-16 ***
mean_diff:symmetric          1   21.8    21.8  438.81 < 2e-16 ***
variance.ratio:equal_size    1    0.0     0.0    0.34  0.5622    
variance.ratio:dists         1    0.0     0.0    0.04  0.8376    
variance.ratio:same_dists    1    0.3     0.3    6.81  0.0092 ** 
variance.ratio:symmetric     1    0.8     0.8   15.25  0.0001 ***
equal_size:dists             1    0.0     0.0    0.38  0.5387    
equal_size:same_dists        1    0.3     0.3    5.38  0.0206 *  
equal_size:symmetric         1    0.3     0.3    6.18  0.0131 *  
dists:same_dists             1    8.8     8.8  177.53 < 2e-16 ***
same_dists:symmetric         1    0.0     0.0    0.13  0.7226    
Residuals                  932   46.2     0.0                    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Going even further, if only two-way interaction with respect to the first factor is of interest (sample_size) then a resolution 4 input can be used, further reducing the total number of conditions to evaluate.

# 7 factors with resolution 4
fr <- FrF2(nfactors=7, resolution=4)

# Create working simulation design given -1/1 combinations
fDesign <- createDesign(sample_size=c(100,200),
                        mean_diff=c(.25, 2),
                        variance.ratio=c(1,4),
                        equal_size=c(TRUE, FALSE),
                        dists=c('norm', 'skew'),
                        same_dists=c(TRUE, FALSE),
                        symmetric=c(TRUE, FALSE),
                        # remove same-normal combo
                        subset = !(symmetric & dists == 'norm'),
                        fractional=fr)
fDesign
# A tibble: 12 × 7
   sample_size mean_diff variance.ratio equal_size dists same_dists symmetric
         <dbl>     <dbl>          <dbl> <lgl>      <chr> <lgl>      <lgl>    
 1         200      0.25              4 FALSE      norm  TRUE       FALSE    
 2         100      0.25              4 TRUE       skew  TRUE       FALSE    
 3         200      2                 1 TRUE       norm  TRUE       FALSE    
 4         100      2                 4 TRUE       norm  FALSE      FALSE    
 5         100      2                 1 FALSE      skew  TRUE       FALSE    
 6         100      0.25              1 FALSE      norm  FALSE      FALSE    
 7         100      0.25              4 FALSE      skew  FALSE      TRUE     
 8         200      2                 4 FALSE      skew  FALSE      FALSE    
 9         200      0.25              1 FALSE      skew  TRUE       TRUE     
10         100      2                 1 TRUE       skew  FALSE      TRUE     
11         200      0.25              1 TRUE       skew  FALSE      FALSE    
12         200      2                 4 TRUE       skew  TRUE       TRUE     
# note the omission of the summarise function
fres <- runSimulation(design=fDesign, replications=20, 
                     generate=Generate, analyse=Analyse, verbose=FALSE)
fres
# A tibble: 240 × 8
   sample_size mean_diff variance.ratio equal_size dists same_dists symmetric
         <dbl>     <dbl>          <dbl> <lgl>      <chr> <lgl>      <lgl>    
 1         200      0.25              4 FALSE      norm  TRUE       FALSE    
 2         200      0.25              4 FALSE      norm  TRUE       FALSE    
 3         200      0.25              4 FALSE      norm  TRUE       FALSE    
 4         200      0.25              4 FALSE      norm  TRUE       FALSE    
 5         200      0.25              4 FALSE      norm  TRUE       FALSE    
 6         200      0.25              4 FALSE      norm  TRUE       FALSE    
 7         200      0.25              4 FALSE      norm  TRUE       FALSE    
 8         200      0.25              4 FALSE      norm  TRUE       FALSE    
 9         200      0.25              4 FALSE      norm  TRUE       FALSE    
10         200      0.25              4 FALSE      norm  TRUE       FALSE    
# ℹ 230 more rows
# ℹ 1 more variable: detect <lgl>
# main effects
mod <- aov(detect ~ ., data=fres)
summary(mod)
                Df Sum Sq Mean Sq F value  Pr(>F)    
sample_size      1   1.35    1.35   20.29 1.1e-05 ***
mean_diff        1   9.60    9.60  144.31 < 2e-16 ***
variance.ratio   1   2.02    2.02   30.32 9.7e-08 ***
equal_size       1   4.80    4.80   72.16 2.4e-15 ***
dists            1   0.17    0.17    2.54  0.1126    
same_dists       1   0.67    0.67   10.15  0.0016 ** 
symmetric        1  13.81   13.81  207.54 < 2e-16 ***
Residuals      232  15.43    0.07                    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# two-way interaction effects for sample_size only
mod2 <- aov(detect ~ (.)^2, data=fres)
summary(mod2)
                            Df Sum Sq Mean Sq F value  Pr(>F)    
sample_size                  1   1.35    1.35   31.73 5.2e-08 ***
mean_diff                    1   9.60    9.60  225.65 < 2e-16 ***
variance.ratio               1   2.02    2.02   47.40 5.5e-11 ***
equal_size                   1   4.80    4.80  112.82 < 2e-16 ***
dists                        1   0.17    0.17    3.97 0.04761 *  
same_dists                   1   0.67    0.67   15.87 9.1e-05 ***
symmetric                    1  13.81   13.81  324.52 < 2e-16 ***
sample_size:mean_diff        1   2.55    2.55   59.99 3.1e-13 ***
sample_size:variance.ratio   1   0.63    0.63   14.69 0.00016 ***
sample_size:equal_size       1   0.06    0.06    1.32 0.25141    
sample_size:dists            1   2.50    2.50   58.76 5.1e-13 ***
Residuals                  228   9.70    0.04                    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

References

  • Montgomery, Douglas C. (2017). Design and Analysis of Experiments (9th ed.). John Wiley & Sons, Inc.