Chapter 48 Fractional Factorials

48.1 FrF2

Package FrF2 has several functions for working with fractioned designs. The basic function is FrF2(). We can use it to construct designs using generators that we specify, or it can select a design for us when we specify number of runs and number of factors.

In the usage here, we want it to find a design for 5 factors in sixteen runs. The output is a data frame with the factor levels (in +1, –1 form). By default, the order of the runs will be randomized.

> library(car)
Loading required package: carData
> library(cfcdae)
Registered S3 method overwritten by 'DoE.base':
  method           from       
  factorize.factor conf.design
> library(FrF2)
Loading required package: DoE.base
Loading required package: grid
Loading required package: conf.design

Attaching package: 'DoE.base'
The following objects are masked from 'package:stats':

    aov, lm
The following object is masked from 'package:graphics':

    plot.design
The following object is masked from 'package:base':

    lengths
> des1 <- FrF2(16,5)
> des1
    A  B  C  D  E
1  -1 -1  1 -1 -1
2  -1  1 -1 -1 -1
3   1 -1 -1  1  1
4  -1  1  1 -1  1
5  -1  1 -1  1  1
6   1  1  1 -1 -1
7   1  1 -1  1 -1
8   1  1  1  1  1
9   1 -1 -1 -1 -1
10 -1 -1 -1 -1  1
11 -1  1  1  1 -1
12  1 -1  1 -1  1
13  1 -1  1  1 -1
14 -1 -1  1  1  1
15 -1 -1 -1  1 -1
16  1  1 -1 -1  1
class=design, type= FrF2 

You can supply factor names if you don’t like A, B, C, etc.

> FrF2(4,3,factor.names=c("Moe","Larry","Curly"))
  Moe Larry Curly
1   1     1     1
2  -1     1    -1
3   1    -1    -1
4  -1    -1     1
class=design, type= FrF2 

You can also name the levels by using a list with named components, each of which contains the level names.

> FrF2(4,3,factor.names=list(speed=c("fast","slow"),temp=c("hot","cold"),
+     time=c("long","short")))
  speed temp  time
1  slow  hot  long
2  slow cold short
3  fast cold  long
4  fast  hot short
class=design, type= FrF2 

In fact, FrF2() has dozens of optional arguments, and we will only scratch the surface of what it can do.

You can get information about a design using the design.info() function. There’s probably more information than you really want.

> design.info(des1)
$type
[1] "FrF2"

$nruns
[1] 16

$nfactors
[1] 5

$factor.names
$factor.names$A
[1] -1  1

$factor.names$B
[1] -1  1

$factor.names$C
[1] -1  1

$factor.names$D
[1] -1  1

$factor.names$E
[1] -1  1


$catlg.name
[1] "catlg"

$catlg.entry
Design:  5-1.1 
   16  runs,  5  factors,  
   Resolution  V 
   Generating columns:  15 
   WLP (3plus):  0 0 1 0 0 ,  10  clear 2fis

$aliased
$aliased$legend
[1] "A=A" "B=B" "C=C" "D=D" "E=E"

$aliased[[2]]
[1] "no aliasing among main effects and 2fis"


$FrF2.version
[1] "2.2-3"

$replications
[1] 1

$repeat.only
[1] FALSE

$randomize
[1] TRUE

$seed
NULL

$creator
FrF2(16, 5)

Most of these designs are going to come from a catalog of designs (like the tables in the back of the book). This command just strips out the catalog information. Some interesting bits include the resolution (here V), and the generating columns.

A \(2^{n-k}\) design is set up as complete factorials in the first \(n-k\) factors. Each of the k additional factors is aliased to an interaction of the first \(n-k\) factors. The “Generating columns” information tells us which columns are the generator; here it is just column 15. These columns are the +1/–1 columns for the various interactions, and they are numbered in standard order beginning with A. Number 15 is the ABCD interaction, so column E (the fifth factor) is generated by ABCD, yielding the overall generator I=ABCDE, giving us resolution V.

> design.info(des1)$catlg.entry
Design:  5-1.1 
   16  runs,  5  factors,  
   Resolution  V 
   Generating columns:  15 
   WLP (3plus):  0 0 1 0 0 ,  10  clear 2fis

Now let’s look at something with more fractioning. Here we have five factors in 8 runs, so this is a quarter fraction that will have two generators.

> des2 <- FrF2(8,5);des2
   A  B  C  D  E
1 -1 -1 -1  1  1
2  1 -1  1 -1  1
3 -1 -1  1  1 -1
4  1  1  1  1  1
5 -1  1  1 -1 -1
6 -1  1 -1 -1  1
7  1 -1 -1 -1 -1
8  1  1 -1  1 -1
class=design, type= FrF2 

FrF2 works by taking a complete factorial in the first \(n-k\) variables and then adding the additional \(k\) variables on using the aliasing. However, it randomizes the order of the runs, so it is often challenging to see the base factorial. If you tell it not to randomize, then you can clearly see the base factorial in standard order.

> FrF2(8,5,randomize=FALSE)
   A  B  C  D  E
1 -1 -1 -1  1  1
2  1 -1 -1 -1 -1
3 -1  1 -1 -1  1
4  1  1 -1  1 -1
5 -1 -1  1  1 -1
6  1 -1  1 -1  1
7 -1  1  1 -1 -1
8  1  1  1  1  1
class=design, type= FrF2 

The generating columns are 3 (AB) and 5 (AC), so the aliasing is based on I = ABD = ACE = BCDE. Note the WLP output. It tells us how many three letter, four letter, and so on words are aliased with I, here 2 three-way interactions and one four-way, as we just saw. It also tells us how many two factor interactions are “clear,” that is, not aliased to main effects or other two factor interactions. In this case, all two factor interactions are aliased to a main effect, so none of them were clear. Just above in des1, all two factor interactions were aliased with three factor interactions, so all (10) of them were clear.

> design.info(des2)$catlg.entry
Design:  5-2.1 
   8  runs,  5  factors,  
   Resolution  III 
   Generating columns:  3 5 
   WLP (3plus):  2 1 0 0 0 ,  0  clear 2fis

Here we ask for a design with generators ACD and BCE, so ABDE is also aliased.

> FrF2(8,generators=c("AC","BC"),randomize=FALSE)
   A  B  C  D  E
1 -1 -1 -1  1  1
2  1 -1 -1 -1  1
3 -1  1 -1  1 -1
4  1  1 -1 -1 -1
5 -1 -1  1 -1 -1
6  1 -1  1  1 -1
7 -1  1  1 -1  1
8  1  1  1  1  1
class=design, type= FrF2.generators 

We can also change signs to get alternate fractions. Here we flip the signs for E.

> FrF2(8,generators=c("AC","-BC"),randomize=FALSE)
   A  B  C  D  E
1 -1 -1 -1  1 -1
2  1 -1 -1 -1 -1
3 -1  1 -1  1  1
4  1  1 -1 -1  1
5 -1 -1  1 -1  1
6  1 -1  1  1  1
7 -1  1  1 -1 -1
8  1  1  1  1 -1
class=design, type= FrF2.generators 

FrF2() always assumes that we have a complete factorial in the first three factors. That means that we could never use ABC as one of the generators. We can get around that by relabeling the factors. In the generator, the AC here really means interaction of first and third factor aliased to the fourth factor; at this stage, FrF2() doesn’t really care how you eventually label the columns. The generators make sure that the fourth column is aliased to the interaction of the first and the third. If the factors are renamed as above, the first, third, and fourth columns are A, B, and C, so we get ABC as one of the aliases of I. Note below that we do not have a full factorial in ABC below; ABC is always positive.

> FrF2(8,generators=c("AC","BC"),factor.nam=c("A","D","C","B","E"),randomize=FALSE)
   A  D  C  B  E
1 -1 -1 -1  1  1
2  1 -1 -1 -1  1
3 -1  1 -1  1 -1
4  1  1 -1 -1 -1
5 -1 -1  1 -1 -1
6  1 -1  1  1 -1
7 -1  1  1 -1  1
8  1  1  1  1  1
class=design, type= FrF2.generators 

We get the same thing using 5 and 6 instead of AC and BC. Recall standard order: A, B, AB, C, AC, BC, ABC, \(\dots\); AC is fifth in that order and BC is sixth. But don’t think about actual factor names, think of it as first column, second column, interaction of first and second, third column, and so on. Then FrF2 adds on the names later.

> FrF2(8,generators=c(5,6),factor.nam=c("A","D","C","B","E"),randomize=FALSE)
   A  D  C  B  E
1 -1 -1 -1  1  1
2  1 -1 -1 -1  1
3 -1  1 -1  1 -1
4  1  1 -1 -1 -1
5 -1 -1  1 -1 -1
6  1 -1  1  1 -1
7 -1  1  1 -1  1
8  1  1  1  1  1
class=design, type= FrF2.generators 

Here is a more interesting design, a \(2^{8-4}\), which is 8 factors in 16 runs for a 1/16 fraction.

> des3 <- FrF2(16,8,randomize=FALSE);des3
    A  B  C  D  E  F  G  H
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
9  -1 -1 -1  1 -1  1  1  1
10  1 -1 -1  1  1 -1 -1  1
11 -1  1 -1  1  1 -1  1 -1
12  1  1 -1  1 -1  1 -1 -1
13 -1 -1  1  1  1  1 -1 -1
14  1 -1  1  1 -1 -1  1 -1
15 -1  1  1  1 -1 -1 -1  1
16  1  1  1  1  1  1  1  1
class=design, type= FrF2 

The generating columns are all three factor interactions [7=ABC, 11=ABD, 13=ACD, 14=BCD], which gives us four factor interactions (and the eight-way interaction) as the aliases of I. Note the resolution IV.

A counts as 1; B counts as 2; C counts as 4; D counts as 8; and so on. Thus ABD is \(1 + 2+ 8 = 11\).

> design.info(des3)$catlg.entry
Design:  8-4.1 
   16  runs,  8  factors,  
   Resolution  IV 
   Generating columns:  7 11 13 14 
   WLP (3plus):  0 14 0 0 0 ,  0  clear 2fis

There are functions for determining the alias structure, but they all use the output of lm(). To use them, we’ll make a dummy response for our design 2 with 8 runs and then fit the model with all main effects and interactions. Clearly, we cannot actually fit all 31 main effects and interactions with only 8 data points.

> des2 <- within(des2,y <- rnorm(8))
> fit2 <- lm(y~A*B*C*D*E,data=des2)

We’ve tried to fit a mean and 31 effects to 8 data points, so only a few things are going to be estimable. The way lm() does things is that it does main effects, then 2fis, etc. It can’t fit AB because it’s aliased with D, it can’t fit AC because it’s aliased with E; it can fit BC, and eventually CD.

> anova(fit2)
Warning in anova.lm(fit2): ANOVA F-tests on an essentially perfect fit are
unreliable
Analysis of Variance Table

Response: y
          Df  Sum Sq Mean Sq F value Pr(>F)
A          1 0.02261 0.02261     NaN    NaN
B          1 1.37322 1.37322     NaN    NaN
C          1 0.06035 0.06035     NaN    NaN
D          1 0.00148 0.00148     NaN    NaN
E          1 0.02042 0.02042     NaN    NaN
B:C        1 0.28338 0.28338     NaN    NaN
C:D        1 0.00076 0.00076     NaN    NaN
Residuals  0 0.00000     NaN               

Lot’s of NAs here.

> fit2

Call:
lm.default(formula = y ~ A * B * C * D * E, data = des2)

Coefficients:
   (Intercept)              A1              B1              C1              D1  
     -0.050267       -0.053159       -0.414309        0.086857        0.013587  
            E1           A1:B1           A1:C1           B1:C1           A1:D1  
     -0.050522              NA              NA       -0.188208              NA  
         B1:D1           C1:D1           A1:E1           B1:E1           C1:E1  
            NA       -0.009716              NA              NA              NA  
         D1:E1        A1:B1:C1        A1:B1:D1        A1:C1:D1        B1:C1:D1  
            NA              NA              NA              NA              NA  
      A1:B1:E1        A1:C1:E1        B1:C1:E1        A1:D1:E1        B1:D1:E1  
            NA              NA              NA              NA              NA  
      C1:D1:E1     A1:B1:C1:D1     A1:B1:C1:E1     A1:B1:D1:E1     A1:C1:D1:E1  
            NA              NA              NA              NA              NA  
   B1:C1:D1:E1  A1:B1:C1:D1:E1  
            NA              NA  

The aliases() function gives us the full table of aliases, except that it does not give the aliases of I. You can find those easily from any row of the output by multiplying by the first element and reducing exponents modulo 2. This gives us I = A:C:E = B:C:D:E = A:B:D.

> aliases(fit2)
                          
 A = C:E = A:B:C:D:E = B:D
 B = C:D:E = A:B:C:E = A:D
 C = A:E = B:D:E = A:B:C:D
 D = B:C:E = A:C:D:E = A:B
 E = B:C:D = A:B:D:E = A:C
 B:C = D:E = A:C:D = A:B:E
 C:D = B:E = A:B:C = A:D:E

The alias() function presents the same information in a different format. There is a column for each fitted term and a row for each aliased term. We then have 1’s wherever an aliased term is aliased to a fitted term.

> alias(fit2)
Model :
y ~ A * B * C * D * E

Complete :
               (Intercept) A1 B1 C1 D1 E1 B1:C1 C1:D1
A1:E1          0           0  0  1  0  0  0     0    
B1:E1          0           0  0  0  0  0  0     1    
C1:E1          0           1  0  0  0  0  0     0    
D1:E1          0           0  0  0  0  0  1     0    
A1:B1:C1       0           0  0  0  0  0  0     1    
A1:B1:D1       1           0  0  0  0  0  0     0    
A1:C1:D1       0           0  0  0  0  0  1     0    
B1:C1:D1       0           0  0  0  0  1  0     0    
A1:B1:E1       0           0  0  0  0  0  1     0    
A1:C1:E1       1           0  0  0  0  0  0     0    
B1:C1:E1       0           0  0  0  1  0  0     0    
A1:D1:E1       0           0  0  0  0  0  0     1    
B1:D1:E1       0           0  0  1  0  0  0     0    
C1:D1:E1       0           0  1  0  0  0  0     0    
A1:B1:C1:D1    0           0  0  1  0  0  0     0    
A1:B1:C1:E1    0           0  1  0  0  0  0     0    
A1:B1:D1:E1    0           0  0  0  0  1  0     0    
A1:C1:D1:E1    0           0  0  0  1  0  0     0    
B1:C1:D1:E1    1           0  0  0  0  0  0     0    
A1:B1:C1:D1:E1 0           1  0  0  0  0  0     0    
A1:B1          0           0  0  0  1  0  0     0    
A1:C1          0           0  0  0  0  1  0     0    
A1:D1          0           0  1  0  0  0  0     0    
B1:D1          0           1  0  0  0  0  0     0    

We can get aliases for bigger designs, but the output can be enormous. (Or, perhaps very tiny if you want to see it all.)

> des3 <- within(des3,y<-rnorm(16))
> fit3 <- lm(y~A*B*C*D*E*F*G*H,data=des3)
> aliases(fit3)
                                                                                                                                                              
 A = B:C:E = B:D:F = C:D:G = E:F:G = D:E:H = C:F:H = B:G:H = A:C:D:E:F = A:B:D:E:G = A:B:C:F:G = A:B:C:D:H = A:B:E:F:H = A:C:E:G:H = A:D:F:G:H = B:C:D:E:F:G:H
 B = A:C:E = A:D:F = D:E:G = C:F:G = C:D:H = E:F:H = A:G:H = B:C:D:E:F = A:B:C:D:G = A:B:E:F:G = A:B:D:E:H = A:B:C:F:H = B:C:E:G:H = B:D:F:G:H = A:C:D:E:F:G:H
 C = A:B:E = D:E:F = A:D:G = B:F:G = B:D:H = A:F:H = E:G:H = A:B:C:D:F = B:C:D:E:G = A:C:E:F:G = A:C:D:E:H = B:C:E:F:H = A:B:C:G:H = C:D:F:G:H = A:B:D:E:F:G:H
 D = A:B:F = C:E:F = A:C:G = B:E:G = B:C:H = A:E:H = F:G:H = A:B:C:D:E = B:C:D:F:G = A:D:E:F:G = A:C:D:F:H = B:D:E:F:H = A:B:D:G:H = C:D:E:G:H = A:B:C:E:F:G:H
 E = A:B:C = C:D:F = B:D:G = A:F:G = A:D:H = B:F:H = C:G:H = A:B:D:E:F = A:C:D:E:G = B:C:E:F:G = B:C:D:E:H = A:C:E:F:H = A:B:E:G:H = D:E:F:G:H = A:B:C:D:F:G:H
 F = A:B:D = C:D:E = B:C:G = A:E:G = A:C:H = B:E:H = D:G:H = A:B:C:E:F = A:C:D:F:G = B:D:E:F:G = B:C:D:F:H = A:D:E:F:H = A:B:F:G:H = C:E:F:G:H = A:B:C:D:E:G:H
 G = A:C:D = B:D:E = B:C:F = A:E:F = A:B:H = C:E:H = D:F:H = A:B:C:E:G = A:B:D:F:G = C:D:E:F:G = B:C:D:G:H = A:D:E:G:H = A:C:F:G:H = B:E:F:G:H = A:B:C:D:E:F:H
 H = B:C:D = A:D:E = A:C:F = B:E:F = A:B:G = C:E:G = D:F:G = A:B:C:E:H = A:B:D:F:H = C:D:E:F:H = A:C:D:G:H = B:D:E:G:H = B:C:F:G:H = A:E:F:G:H = A:B:C:D:E:F:G
 A:B = D:F = G:H = B:C:D:G = A:D:E:G = A:C:F:G = B:E:F:G = A:C:D:H = B:D:E:H = B:C:F:H = A:E:F:H = A:B:C:D:E:F = A:B:C:E:G:H = A:B:D:F:G:H = C:D:E:F:G:H = C:E
 A:C = D:G = F:H = B:C:D:F = A:D:E:F = A:B:F:G = C:E:F:G = A:B:D:H = C:D:E:H = B:C:G:H = A:E:G:H = A:B:C:D:E:G = A:B:C:E:F:H = A:C:D:F:G:H = B:D:E:F:G:H = B:E
 B:C = F:G = D:H = A:C:D:F = B:D:E:F = A:B:D:G = C:D:E:G = A:B:F:H = C:E:F:H = A:C:G:H = B:E:G:H = A:B:C:E:F:G = A:B:C:D:E:H = B:C:D:F:G:H = A:D:E:F:G:H = A:E
 A:D = B:F = C:G = E:H = B:C:D:E = A:C:E:F = A:B:E:G = D:E:F:G = A:B:C:H = C:D:F:H = B:D:G:H = A:F:G:H = A:B:C:D:F:G = A:B:D:E:F:H = A:C:D:E:G:H = B:C:E:F:G:H
 B:D = A:F = E:G = C:H = A:C:D:E = B:C:E:F = A:B:C:G = C:D:F:G = A:B:E:H = D:E:F:H = A:D:G:H = B:F:G:H = A:B:D:E:F:G = A:B:C:D:F:H = B:C:D:E:G:H = A:C:E:F:G:H
 C:D = E:F = A:G = B:H = A:B:D:E = A:B:C:F = B:C:E:G = B:D:F:G = A:C:E:H = A:D:F:H = D:E:G:H = C:F:G:H = A:C:D:E:F:G = B:C:D:E:F:H = A:B:C:D:G:H = A:B:E:F:G:H
 D:E = C:F = B:G = A:H = A:B:C:D = A:B:E:F = A:C:E:G = A:D:F:G = B:C:E:H = B:D:F:H = C:D:G:H = E:F:G:H = B:C:D:E:F:G = A:C:D:E:F:H = A:B:D:E:G:H = A:B:C:F:G:H

Missing data change the aliasing structure, and usually make it much more complicated. Here we will just make a second response that has one missing value.

> des2 <- within(des2,{yb <- y;yb[1]<-NA})
> fit2b <- lm(yb~A*B*C*D*E,data=des2)
> alias(fit2b)
Model :
yb ~ A * B * C * D * E

Complete :
               (Intercept) A1 B1 C1 D1 E1 B1:C1
A1:D1           0           0  1  0  0  0  0   
B1:D1           0           1  0  0  0  0  0   
C1:D1           1          -1 -1 -1  1  1  1   
A1:E1           0           0  0  1  0  0  0   
B1:E1           1          -1 -1 -1  1  1  1   
C1:E1           0           1  0  0  0  0  0   
D1:E1           0           0  0  0  0  0  1   
A1:B1:C1        1          -1 -1 -1  1  1  1   
A1:B1:D1        1           0  0  0  0  0  0   
A1:C1:D1        0           0  0  0  0  0  1   
B1:C1:D1        0           0  0  0  0  1  0   
A1:B1:E1        0           0  0  0  0  0  1   
A1:C1:E1        1           0  0  0  0  0  0   
B1:C1:E1        0           0  0  0  1  0  0   
A1:D1:E1        1          -1 -1 -1  1  1  1   
B1:D1:E1        0           0  0  1  0  0  0   
C1:D1:E1        0           0  1  0  0  0  0   
A1:B1:C1:D1     0           0  0  1  0  0  0   
A1:B1:C1:E1     0           0  1  0  0  0  0   
A1:B1:D1:E1     0           0  0  0  0  1  0   
A1:C1:D1:E1     0           0  0  0  1  0  0   
B1:C1:D1:E1     1           0  0  0  0  0  0   
A1:B1:C1:D1:E1  0           1  0  0  0  0  0   
A1:B1           0           0  0  0  1  0  0   
A1:C1           0           0  0  0  0  1  0   

Speaking of messy, before moving on to data, let’s look at another class of fractional factorial designs called Plackett-Burman designs. They are resolution III designs in k=N-1 factors in N data points when N is a multiple of 4. In this example, we look at 11 factors in 12 runs.

In principle we could include all 4047 main effects and interactions, but I just went up to 2 factor interactions, as they are sufficient to illustrate the point.

> A <- factor(c(2,2,1,2,2,2,1,1,1,2,1,1))
> B <- factor(c(1,2,2,1,2,2,2,1,1,1,2,1))
> C <- factor(c(2,1,2,2,1,2,2,2,1,1,1,1))
> D <- factor(c(1,2,1,2,2,1,2,2,2,1,1,1))
> E <- factor(c(1,1,2,1,2,2,1,2,2,2,1,1))
> F <- factor(c(1,1,1,2,1,2,2,1,2,2,2,1))
> G <- factor(c(2,1,1,1,2,1,2,2,1,2,2,1))
> H <- factor(c(2,2,1,1,1,2,1,2,2,1,2,1))
> I <- factor(c(2,2,2,1,1,1,2,1,2,2,1,1))
> J <- factor(c(1,2,2,2,1,1,1,2,1,2,2,1))
> K <- factor(c(2,1,2,2,2,1,1,1,2,1,2,1))
> pb <- data.frame(A,B,C,D,E,F,G,H,I,J,K,y=rnorm(12))
> outpb <- lm(y~(A+B+C+D+E+F+G+H+I+J+K)^2,data=pb)

Here we look at aliasing. Every second order interaction is aliased to all main effects not occurring in the interaction. Thus you had better be really, really sure about not having any sizable interactions when you use a PB design, because a single non-ignorable interaction screws up all of your main effects estimates, not just one.

> alias(outpb)
Model :
y ~ (A + B + C + D + E + F + G + H + I + J + K)^2

Complete :
      (Intercept) A1   B1   C1   D1   E1   F1   G1   H1   I1   J1   K1  
A1:B1    0           0    0  1/3 -1/3 -1/3  1/3  1/3 -1/3  1/3  1/3  1/3
A1:C1    0           0  1/3    0  1/3  1/3 -1/3  1/3 -1/3  1/3  1/3 -1/3
A1:D1    0           0 -1/3  1/3    0  1/3  1/3  1/3  1/3  1/3 -1/3 -1/3
A1:E1    0           0 -1/3  1/3  1/3    0 -1/3 -1/3  1/3  1/3  1/3  1/3
A1:F1    0           0  1/3 -1/3  1/3 -1/3    0  1/3  1/3  1/3 -1/3  1/3
A1:G1    0           0  1/3  1/3  1/3 -1/3  1/3    0  1/3 -1/3  1/3 -1/3
A1:H1    0           0 -1/3 -1/3  1/3  1/3  1/3  1/3    0 -1/3  1/3  1/3
A1:I1    0           0  1/3  1/3  1/3  1/3  1/3 -1/3 -1/3    0 -1/3  1/3
A1:J1    0           0  1/3  1/3 -1/3  1/3 -1/3  1/3  1/3 -1/3    0  1/3
A1:K1    0           0  1/3 -1/3 -1/3  1/3  1/3 -1/3  1/3  1/3  1/3    0
B1:C1    0         1/3    0    0  1/3 -1/3 -1/3  1/3  1/3 -1/3  1/3  1/3
B1:D1    0        -1/3    0  1/3    0  1/3  1/3 -1/3  1/3 -1/3  1/3  1/3
B1:E1    0        -1/3    0 -1/3  1/3    0  1/3  1/3  1/3  1/3  1/3 -1/3
B1:F1    0         1/3    0 -1/3  1/3  1/3    0 -1/3 -1/3  1/3  1/3  1/3
B1:G1    0         1/3    0  1/3 -1/3  1/3 -1/3    0  1/3  1/3  1/3 -1/3
B1:H1    0        -1/3    0  1/3  1/3  1/3 -1/3  1/3    0  1/3 -1/3  1/3
B1:I1    0         1/3    0 -1/3 -1/3  1/3  1/3  1/3  1/3    0 -1/3  1/3
B1:J1    0         1/3    0  1/3  1/3  1/3  1/3  1/3 -1/3 -1/3    0 -1/3
B1:K1    0         1/3    0  1/3  1/3 -1/3  1/3 -1/3  1/3  1/3 -1/3    0
C1:D1    0         1/3  1/3    0    0  1/3 -1/3 -1/3  1/3  1/3 -1/3  1/3
C1:E1    0         1/3 -1/3    0  1/3    0  1/3  1/3 -1/3  1/3 -1/3  1/3
C1:F1    0        -1/3 -1/3    0 -1/3  1/3    0  1/3  1/3  1/3  1/3  1/3
C1:G1    0         1/3  1/3    0 -1/3  1/3  1/3    0 -1/3 -1/3  1/3  1/3
C1:H1    0        -1/3  1/3    0  1/3 -1/3  1/3 -1/3    0  1/3  1/3  1/3
C1:I1    0         1/3 -1/3    0  1/3  1/3  1/3 -1/3  1/3    0  1/3 -1/3
C1:J1    0         1/3  1/3    0 -1/3 -1/3  1/3  1/3  1/3  1/3    0 -1/3
C1:K1    0        -1/3  1/3    0  1/3  1/3  1/3  1/3  1/3 -1/3 -1/3    0
D1:E1    0         1/3  1/3  1/3    0    0  1/3 -1/3 -1/3  1/3  1/3 -1/3
D1:F1    0         1/3  1/3 -1/3    0  1/3    0  1/3  1/3 -1/3  1/3 -1/3
D1:G1    0         1/3 -1/3 -1/3    0 -1/3  1/3    0  1/3  1/3  1/3  1/3
D1:H1    0         1/3  1/3  1/3    0 -1/3  1/3  1/3    0 -1/3 -1/3  1/3
D1:I1    0         1/3 -1/3  1/3    0  1/3 -1/3  1/3 -1/3    0  1/3  1/3
D1:J1    0        -1/3  1/3 -1/3    0  1/3  1/3  1/3 -1/3  1/3    0  1/3
D1:K1    0        -1/3  1/3  1/3    0 -1/3 -1/3  1/3  1/3  1/3  1/3    0
E1:F1    0        -1/3  1/3  1/3  1/3    0    0  1/3 -1/3 -1/3  1/3  1/3
E1:G1    0        -1/3  1/3  1/3 -1/3    0  1/3    0  1/3  1/3 -1/3  1/3
E1:H1    0         1/3  1/3 -1/3 -1/3    0 -1/3  1/3    0  1/3  1/3  1/3
E1:I1    0         1/3  1/3  1/3  1/3    0 -1/3  1/3  1/3    0 -1/3 -1/3
E1:J1    0         1/3  1/3 -1/3  1/3    0  1/3 -1/3  1/3 -1/3    0  1/3
E1:K1    0         1/3 -1/3  1/3 -1/3    0  1/3  1/3  1/3 -1/3  1/3    0
F1:G1    0         1/3 -1/3  1/3  1/3  1/3    0    0  1/3 -1/3 -1/3  1/3
F1:H1    0         1/3 -1/3  1/3  1/3 -1/3    0  1/3    0  1/3  1/3 -1/3
F1:I1    0         1/3  1/3  1/3 -1/3 -1/3    0 -1/3  1/3    0  1/3  1/3
F1:J1    0        -1/3  1/3  1/3  1/3  1/3    0 -1/3  1/3  1/3    0 -1/3
F1:K1    0         1/3  1/3  1/3 -1/3  1/3    0  1/3 -1/3  1/3 -1/3    0
G1:H1    0         1/3  1/3 -1/3  1/3  1/3  1/3    0    0  1/3 -1/3 -1/3
G1:I1    0        -1/3  1/3 -1/3  1/3  1/3 -1/3    0  1/3    0  1/3  1/3
G1:J1    0         1/3  1/3  1/3  1/3 -1/3 -1/3    0 -1/3  1/3    0  1/3
G1:K1    0        -1/3 -1/3  1/3  1/3  1/3  1/3    0 -1/3  1/3  1/3    0
H1:I1    0        -1/3  1/3  1/3 -1/3  1/3  1/3  1/3    0    0  1/3 -1/3
H1:J1    0         1/3 -1/3  1/3 -1/3  1/3  1/3 -1/3    0  1/3    0  1/3
H1:K1    0         1/3  1/3  1/3  1/3  1/3 -1/3 -1/3    0 -1/3  1/3    0
I1:J1    0        -1/3 -1/3  1/3  1/3 -1/3  1/3  1/3  1/3    0    0  1/3
I1:K1    0         1/3  1/3 -1/3  1/3 -1/3  1/3  1/3 -1/3    0  1/3    0
J1:K1    0         1/3 -1/3 -1/3  1/3  1/3 -1/3  1/3  1/3  1/3    0    0

48.2 Analysis

Chen(2003) via Lye (2019) has a \(2^{6-2}\) design on the physical cleaning of reverse osmosis membranes. The six factors are:

  • Production interval between physical cleaning (.5 or 3 hours)
  • Duration of forward flush (1 or 5 minutes)
  • Duration of backwash (1 or 5 minutes)
  • Pressure during forward flush (172 or 345 kPa)
  • Type of water used (RO permeate or Tap water)
  • Sequence (backward then forward, or forward then backward)

The response of interest is wash water usage (WWU).

Here we read in the data and determine that the aliasing is
I = ACDF = ABEF = BCDE.

> washwater <- read.csv("WashWater.csv",header=TRUE)
> washwater <- within(washwater,{A<-factor(A);B<-factor(B);C<-factor(C);D<-factor(D);E<-factor(E);F<-factor(F)})
> fit.wwu <- lm(WWU~A*B*C*D*E*F,data=washwater)
> aliases(fit.wwu)
                              
 A = C:D:F = B:E:F = A:B:C:D:E
 B = C:D:E = A:E:F = A:B:C:D:F
 C = B:D:E = A:D:F = A:B:C:E:F
 D = B:C:E = A:C:F = A:B:D:E:F
 E = B:C:D = A:B:F = A:C:D:E:F
 F = A:C:D = A:B:E = B:C:D:E:F
 A:B = A:C:D:E = B:C:D:F = E:F
 A:C = A:B:D:E = B:C:E:F = D:F
 B:C = A:B:D:F = A:C:E:F = D:E
 A:D = A:B:C:E = B:D:E:F = C:F
 B:D = A:B:C:F = A:D:E:F = C:E
 C:D = A:B:C:D:E:F = B:E = A:F
 A:E = A:B:C:D = C:D:E:F = B:F
 A:B:C = A:D:E = B:D:F = C:E:F
 A:B:D = A:C:E = B:C:F = D:E:F

To do the analysis, we fit to a full factorial. There are a bunch available, but let’s try the last four variables.

> fit2.wwu <- lm(WWU~C*D*E*F,data=washwater)
> TwoSeriesPlots(fit2.wwu)
B set to 3000
Two series plot for wash water usage data

Figure 48.1: Two series plot for wash water usage data

The big effects are CDF, C, CDE, DF, and EF. That’s a bizarre set, but this is a fraction, so check the aliasing. If you go through the aliasing table above, you find that CDF is aliased to A; CDE is aliased to B; DF is aliased to AC; and EF is aliased to AB. So what we are really looking at is A, B, C, AB and AC, which makes more sense.

Let’s refit with just the big effects and take a look at the residuals. We see nonconstant variance.

> fit3.wwu <- lm(WWU~A*B*C,data=washwater)
> plot(fit3.wwu,which=1)
Resdual plot of wash water data reduced model

Figure 48.2: Resdual plot of wash water data reduced model

Box-Cox says try a log. Note that this BC plot is unusually asymmetric.

> boxCox(fit3.wwu)
Box-Cox plot of wash water data reduced model

Figure 48.3: Box-Cox plot of wash water data reduced model

Transform, refit, and look at the plot.

The AB and AC interactions no longer look large; that is not too surprising. The BC interaction (aliased to DE) looking large (on the log scale) was, I have to confess, not something I was expecting.

> lWWU <- log(washwater$WWU)
> fit4.wwu <- lm(lWWU~C*D*E*F,data=washwater)
> TwoSeriesPlots(fit4.wwu)
B set to 3000
Two series plot for logged wash water data

Figure 48.4: Two series plot for logged wash water data

Fit the reduced model and look at the residuals; they’re somewhat better.

> fit5.wwu <- lm(lWWU~A+B+C+B:C,data=washwater)
> plot(fit5.wwu,which=1)
Residual plot for reduced model of log wash water

Figure 48.5: Residual plot for reduced model of log wash water

Now let’s look at that interaction. The standard errors are perhaps not totally trustworthy (we’ve put the small stuff into error, so they may be somewhat too short), but the interaction is clear and clearly bigger than noise: B has a much larger effect at the low level of C (and vice versa).

> anova(fit5.wwu)
Analysis of Variance Table

Response: lWWU
          Df Sum Sq Mean Sq F value    Pr(>F)    
A          1 8.2247  8.2247 709.173 2.432e-11 ***
B          1 2.0049  2.0049 172.869 4.524e-08 ***
C          1 1.9248  1.9248 165.961 5.591e-08 ***
B:C        1 0.3776  0.3776  32.561 0.0001369 ***
Residuals 11 0.1276  0.0116                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> with(washwater,interactplot(B,C,lWWU,sigma2=.0116,df=11,confidence = .95))
B:C interaction plot for wash water data

Figure 48.6: B:C interaction plot for wash water data