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

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)

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)

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

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)

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))

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