Resolution 5 Fractional Factorial Designs

The Problem

At work, I needed to create a resolution V (R5) fractional factorial design for use in a study of the relative effectiveness of various modernization programs. The computer network on which I was working was an isolated network completely disconnected from the internet with specific versions of software installed to support the type of work done. When I’ve needed to create an R5 design in the past, I either used the FrF2 R package or a Ruby script written by Dr. Paul Sanchez at the Naval Postgraduate School. Neither of those were viable options on the isolated computer network, so I needed to create an R5 design for 14 factors from scratch. Since base R was available for this task, that seemed the path of least resistance.

Resolution V Designs

So why did I need to create an R5 design, and what is one? The study involved 14 modernization programs, and I wanted to evaluate the contribution of each modernization program compared to its current equivalent. That means I had 14 factors with two levels each: current and future. The approach was to represent each of the factors at both levels in a computer simulation, query the simulation results for one or more measure of effectiveness, and fit one or more machine learning models to the results to understand the relationship between the factors and the measures of effectiveness.

One approach to generate these results is to run the simulation for each unique combination of factors and their levels. 14 factors at 2 levels each would require simulation runs. A benefit of this design (called a full factorial design) is that it allows for the evaluation of main effects, 2-way interactions, 3-way interactions, and all the way up to 13-way interactions. Full factorial designs are easy to create in R using the expand.grid() function. An example for thee factors is shown below, and this is often referred to as a design matrix. Note the use of -1 and 1 to represent the two factor levels instead of 0 and 1 will be explained later. When evaluating the properties of the design itself, -1 and 1 is necessary (and will be demonstrated later), but we re-code to 0 and 1 when applying machine learning methods to the results.

expand.grid(
  f1 = c(-1, 1),
  f2 = c(-1, 1),
  f3 = c(-1, 1)
  )

##   f1 f2 f3
## 1 -1 -1 -1
## 2  1 -1 -1
## 3 -1  1 -1
## 4  1  1 -1
## 5 -1 -1  1
## 6  1 -1  1
## 7 -1  1  1
## 8  1  1  1

For this study, each simulation took approximately one hour of computing time, and the network was fairly small, so there were a limited number of machines on which to farm out all of these runs. Each simulation also produced approximately 10Gb of output. Unfortunately, the computer network simply didn’t have the computing power or storage space to handle all 16,384 runs.

Clearly, we needed to reduce the number of simulation runs to something feasible. That’s what an R5 design does for us, but it comes at a cost. If we were willing to give up the ability to evaluate 3-way and higher order interactions (we were), then we could drastically reduce the number of simulation runs. For example, to evaluate main effects and all 2-way interactions for 14 factors, we know we need at least simulation runs to provide the necessary degrees of freedom - a huge decrease from 16,384. We might need more than 106 runs, but not orders of magnitude more.

So we know we needed a design matrix with at least 106 rows, but how do we create it? Montgomery (2013) discusses a method that uses the smallest full factorial design with more rows than needed for the R5 design as a basis. The smallest full factorial design with at least 106 rows is , so we have a design matrix of size 128 x 128. From this matrix, we would then select 14 columns to represent our 14 factors. Montgomery (2013) specifies which columns to select for R5 designs but only up to 10 factors. One option for creating our design would be to continue with this approach of using full factorial designs as a basis. Dr. Sanchez’s Ruby script uses a different option, which piqued my interest. Instead of a full factorial design as the basis, his script uses a naturally ordered Walsh matrix (also referred to as a Hadamard matrix) as the basis.

Walsh Matrices

An internet search will produce plenty of options to learn all about Walsh matrices, but I’ll boil that down into their properties that matter for this purpose.

• As with full factorial designs, they are square with dimensions of .

• They consist entirely of binary values.

• All columns (and less importantly, rows) are orthogonal.

• They are easy to generate.

The smallest Walsh matrix is 2x2 with the following values.

w <- matrix(c(1, 1, 1, -1), nrow = 2)
w

##      [,1] [,2]
## [1,]    1    1
## [2,]    1   -1

The Kronecker product of two 2x2 Walsh matrices creates a 4x4 Walsh matrix.

big_w <- w %x% w
big_w

##      [,1] [,2] [,3] [,4]
## [1,]    1    1    1    1
## [2,]    1   -1    1   -1
## [3,]    1    1   -1   -1
## [4,]    1   -1   -1    1

If you’re like me and have never heard of the Kronecker product before, I’ve highlighted and divided the matrix into sections to illustrate what it does. Take the upper left value in the right hand matrix and multiply it by the entire left hand matrix. Place that matrix in the upper left quadrant of the new matrix big_w. Then cycle through each value in the right hand matrix in the same manner.

V1	V2	V3	V4
1	1	1	1
1	-1	1	-1
1	1	-1	-1
1	-1	-1	1

To create an 8x8 matrix, take the Kronecker product of the 4x4 and 2x2 Walsh matrices. Repeating this process will generate increasingly large matrices bounded only by available RAM.

big_w <- big_w %x% w
big_w

##      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## [1,]    1    1    1    1    1    1    1    1
## [2,]    1   -1    1   -1    1   -1    1   -1
## [3,]    1    1   -1   -1    1    1   -1   -1
## [4,]    1   -1   -1    1    1   -1   -1    1
## [5,]    1    1    1    1   -1   -1   -1   -1
## [6,]    1   -1    1   -1   -1    1   -1    1
## [7,]    1    1   -1   -1   -1   -1    1    1
## [8,]    1   -1   -1    1   -1    1    1   -1

Determine Valid Matrix Size and Column Indices

Let’s start with the smallest case possible: create an R5 design for two factors. The first step is to identify the theoretical minimum number of runs using the method described earlier: intercept + # main effects + # interactions. Then start with the smallest Walsh matrix with at least that many rows. For two factors, that’s 4 rows, so we need a 4x4 Walsh matrix. We also know we’re going to drop the first column, so we have 3 columns to work with for our 2 factors.

# minimum number of rows
1 + 2 + choose(2, 2)

## [1] 4

# get Walsh matrix and display only last three columns
big_w <- w %x% w
big_w

##      [,1] [,2] [,3] [,4]
## [1,]    1    1    1    1
## [2,]    1   -1    1   -1
## [3,]    1    1   -1   -1
## [4,]    1   -1   -1    1

We don’t need to check if main effects are aliased because we know they aren’t based on the properties of Walsh matrices. We just need to check whether main effects are aliased with 2-way interactions. Later we’ll also need to check whether 2-way interactions are aliased with any other 2-way interactions. We can skip that for 2 factors since there’s only one 2-way interaction. So, we have three possible candidates for our 2 factors: columns 2, 3, and 4. Let’s just start at the left and work to the right as needed. We’ll select columns 2 and 3 and check the correlation between the two main effects with the interaction term. If no correlation exists, then we’ll keep track of the size of the matrix we tried and the indices of the columns that produced no correlation.

# get the interaction between the candidate columns
int_term <- big_w[, 2] * big_w[, 3]

# check for correlation of both main effects and the interaction term
cor(big_w[, 2], int_term)

## [1] 0

cor(big_w[, 3], int_term)

## [1] 0

Since there is no correlation, we know that for 2 factors, we can use columns 2 and 3 from a 4x4 matrix. If we weren’t able to find a combination of columns that worked, we’d increase the size of the underlying Walsh matrix and repeat the process. The pwr vector will keep track of the number of times we needs to use the Kronecker product to generate the correct sized matrix. For 2 factors and a 4x4 matrix, we use the kronecker product once, so we add a 1 to the vector at an index of 2.

idx <- c(2, 3)
pwr <- c(0, 1)

Adding Factors

Now that we know the indices and matrix size for two factors, we can use this same algorithm to successively add one more factor at a time.

# minimum number of rows
1 + 3 + choose(3, 2)

## [1] 7

So we need at least an 8x8 matrix.

big_w <- big_w %x% w
dim(big_w)

## [1] 8 8

We’re already using columns 2 and 3, so for our third factor, we’ll start with column 4 as a candidate and check for non-zero correlation. If it fails the correlation checks, we’ll move on to column 5, and continue until we find a valid column. If we find a valid column, we’ll record the column index number.

fn <- function(x){m[, x[1]] * m[, x[2]]}

# loop through the candidate columns
for (i in (1+max(idx)):ncol(big_w)){
  # create a new dataframe with just the columns of interest
  m <- big_w[, c(idx, i)]
  # check for correlation. If none found, save the column number and stop the loop.
  if (all(colSums(cor(cbind(m[, ncol(m)], combn(1:(length(idx)+1), 2, fn)))) == 1)){
    idx <- c(idx, i)
    break
  }
}
idx

## [1] 2 3 5

Looks like column 4 failed the test but column 5 passed. We found a valid column, so we’ll also record the number of times we invoked the kronecker product (2) in the pwr vector power at index 3 (because it corresponds with 3 factors).

pwr <- c(pwr, 2)

Continuing On

From here forward, it’s the same process applied to each new column. However, once the number of indices in the idx vector gets into the 30s, the algorithm as written starts slowing down drastically. By avoiding checking correlations that have been checked during earlier iterations and improving the efficiency of my code, I was able to generate R5 matrices for up to 70 factors. At that point, the base matrix dimensions were approximately 32,000 by 32,000 which required 8GB of RAM. My personal computer has 16GB of RAM, and so due to R’s copy on modify behavior, I was unable to create a matrix of that size. I’ve yet to be involved in a DOE study where more than 70 factors were being considered, so I’ll just accept that limitation for now. The full algorithm for up to 20 factors is shown below.

# function to create a Walsh matrix of a certain size.
# Since I iterate from 2:p, I get the size of the Walsh matrix as 2^p
res5 <- function(fct, p){
  a <- matrix(c(1,1,1,-1), nrow=2)
  w <- matrix(c(1,1,1,-1), nrow=2)
  for (i in 2:p){w <- w %x% a}
  w
}

idx <- c(2, 3)
pwr <- c(0, 1)

# function to check aliasing between main effects and two-way interactions
one_two <- function(i) {sum(m[, ncol(m)] == m[, int1[1, i]] * m[, int1[2, i]]) == nrow(m) |
    sum(m[, ncol(m)] == -(m[, int1[1, i]] * m[, int1[2, i]])) == nrow(m)}

# iterate through new factors
for (f in 3:20){
  # get the minimum size of the Walsh matrix
  new_p <- max(ceiling(log2(1 + f + choose(f, 2))), pwr[length(pwr)])
  # try the smallest matrix size first, then try one bigger
  for (ps in c(new_p, new_p+1)){
    p_worked <- FALSE
    dm <- res5(f, ps)
    # get the various combinations of interactions
    int1 <- combn(1:length(idx), 2)
    int2 <- combn(1:(length(idx)+1), 2)
    # filter out the interactions we've already checked
    int_n <- int2[, int2[2, ]==f]
    # iterate through the new potential matrix columns
    for (k in (max(idx)+1):ncol(dm)){
      m <- dm[, c(idx, k)]
      # check new factor against 2-way interaction terms
      keepchecking <- TRUE
      if(any(1:ncol(int1) %>% purrr::map_lgl(function(x) one_two(x)))) {next}
      # check 2-way interactions vs. 2-way interactions
      if (keepchecking){
        for (j in 1:ncol(int1)){
          for (i in 1:ncol(int_n)){
            if(cor(m[, int1[1, j]] * m[, int1[2, j]], m[, int_n[1, i]] * m[, int_n[2, i]]) > 0){
              keepchecking <- FALSE
              break
            }
          }
          if (!keepchecking){break}
        }
      }
      if (keepchecking){
        idx <- c(idx, k)
        pwr <- c(pwr, ps)
        p_worked <- TRUE
        break}
    }
    if (p_worked){break}
  }
}

The column indices are:

idx

##  [1]   2   3   5   9  16  17  33  52  65  86 107 129 151 172 220 238 248 257 280
## [20] 298

And the associated sizes of the Walsh matrices are:

pwr

##  [1] 0 1 3 4 4 5 6 6 7 7 7 8 8 8 8 8 8 9 9 9

Now we have everything needed to create R5 designs for up to 20 factors. We’ll wrap everything up in a new function and replace all -1’s with 0’s for convenience. Then we’ll display the first few rows of the design matrix for 14 factors.

r5 <- function(fct){
  idx <- c(2, 3, 5, 9, 16, 17, 33, 52, 65, 86, 107, 129, 151, 172, 220,
           238, 248, 257, 280, 298)
  pwr <- c(0, 1, 3, 4, 4, 5, 6, 6, 7, 7, 7, 8, 8, 8, 8, 8, 8, 9, 9, 9)

  a <- matrix(c(1,1,1,-1), nrow=2)
  w <- matrix(c(1,1,1,-1), nrow=2)
  for (i in 2:pwr[fct]){w <- w %x% a}
  # select only the columns we need
  w <- w[, idx[1:fct]]
  # replace -1 with 0
  w[w == -1] <- 0
  w
}

head(r5(14))

##      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14]
## [1,]    1    1    1    1    1    1    1    1    1     1     1     1     1     1
## [2,]    0    1    1    1    0    1    1    0    1     0     1     1     1     0
## [3,]    1    0    1    1    0    1    1    0    1     1     0     1     0     0
## [4,]    0    0    1    1    1    1    1    1    1     0     0     1     0     1
## [5,]    1    1    0    1    0    1    1    1    1     0     1     1     0     1
## [6,]    0    1    0    1    1    1    1    0    1     1     1     1     0     0

The number of runs required for this design is simply nrow(r5(14)) = 256.

I am currently working with a colleague to develop an R package that uses the methodology described above to generate Resolution III, V, and VII fractional factorial designs. The package will also provide the ability to generate nearly orthogonal Latin hypercube designs.

Reference

Montgomery, Douglas C. 2013. Design and Analysis of Experiments. John Wiley & Sons.