R Tools for Two-Level Designs
Working with two-level factorial designs is easier with the right tools.
Basic design
This simple R function:
twoLevel = function(k, names = LETTERS[LETTERS != "I"][1:k]) {
x = NULL;
for (j in 1:k)
x = rbind(cbind(x, -1), cbind(x, +1));
dimnames(x)[[2]] = names;
as.data.frame(x);
}
creates a basic 2k design in the standard order; for example 23:
> twoLevel(3); A B C 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
Note that the runs would need to be randomized to produce a final design.
A data-frame suitable for analysis using lm()
or aov() may be prepared by adding the
response:
> filtration = twoLevel(4);
> filtration$Rate = c(45, 71, 48, 65, 68, 60, 80, 65,
43, 100, 45, 104, 75, 86, 70, 96);
> filtration
A B C D Rate
1 -1 -1 -1 -1 45
2 1 -1 -1 -1 71
3 -1 1 -1 -1 48
4 1 1 -1 -1 65
5 -1 -1 1 -1 68
6 1 -1 1 -1 60
7 -1 1 1 -1 80
8 1 1 1 -1 65
9 -1 -1 -1 1 43
10 1 -1 -1 1 100
11 -1 1 -1 1 45
12 1 1 -1 1 104
13 -1 -1 1 1 75
14 1 -1 1 1 86
15 -1 1 1 1 70
16 1 1 1 1 96
Blocking and Replication
Call twoLevel() multiple times to create a
blocked full factorial design; for example, 22 in
three blocks:
> block1 = block2 = block3 = twoLevel(2);
> block1$Rep = "Rep1";
> block2$Rep = "Rep2";
> block3$Rep = "Rep3";
> rbind(block1, block2, block3);
A B Rep
1 -1 -1 Rep1
2 1 -1 Rep1
3 -1 1 Rep1
4 1 1 Rep1
5 -1 -1 Rep2
6 1 -1 Rep2
7 -1 1 Rep2
8 1 1 Rep2
9 -1 -1 Rep3
10 1 -1 Rep3
11 -1 1 Rep3
12 1 1 Rep3
If the design is really replicated, not in blocks, presenting the rows in standard order may be convenient:
> a = rbind(block1, block2, block3);
> a[order(a$B, a$A),]
A B Rep
1 -1 -1 Rep1
5 -1 -1 Rep2
9 -1 -1 Rep3
2 1 -1 Rep1
6 1 -1 Rep2
10 1 -1 Rep3
3 -1 1 Rep1
7 -1 1 Rep2
11 -1 1 Rep3
4 1 1 Rep1
8 1 1 Rep2
12 1 1 Rep3
Blocking with Confounding
If blocks contain too few runs for the full basic design, one or more effects will be confounded with blocks; for example, 24 in 2 blocks each of 8 runs, with ABCD confounded with blocks:
> a = twoLevel(4);
> a$Block = a$A * a$B * a$C * a$D;
> a
A B C D Block
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
The design would usually be shown with the runs sorted into blocks:
> a[order(a$Block, a$C, a$B, a$A),]
A B C D Block
9 -1 -1 -1 1 -1
2 1 -1 -1 -1 -1
3 -1 1 -1 -1 -1
12 1 1 -1 1 -1
5 -1 -1 1 -1 -1
14 1 -1 1 1 -1
15 -1 1 1 1 -1
8 1 1 1 -1 -1
1 -1 -1 -1 -1 1
10 1 -1 -1 1 1
11 -1 1 -1 1 1
4 1 1 -1 -1 1
13 -1 -1 1 1 1
6 1 -1 1 -1 1
7 -1 1 1 -1 1
16 1 1 1 1 1
Use the alias() function to detect confounding;
the function requires a response in the formula, but any
values, such as all zeroes, will do:
> a$y = 0;
> alias(y ~ Block + A * B * C * D, a)
Model :
y ~ Block + A * B * C * D
Complete :
(Intercept) Block A B C D A:B A:C B:C A:D B:D C:D A:B:C A:B:D A:C:D
A:B:C:D 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0
B:C:D
A:B:C:D 0
The alias table has only a single row (but wrapped), and it
confirms that the A:B:C:D interaction is
completely confounded with blocks.
To use more, smaller blocks, more effects will be confounded with blocks; construct unique blocks from chosen confounded effects; for example, 24 in 4 blocks each of 4 runs, with ABC and ABD (and hence also CD) confounded with blocks:
> a$Block = as.factor(paste((1 + a$A * a$B * a$C) / 2,
(1 + a$A * a$B * a$D) / 2, sep = ""));
> a[order(a$Block),]
A B C D Block
1 -1 -1 -1 -1 00
4 1 1 -1 -1 00
14 1 -1 1 1 00
15 -1 1 1 1 00
6 1 -1 1 -1 01
7 -1 1 1 -1 01
9 -1 -1 -1 1 01
12 1 1 -1 1 01
5 -1 -1 1 -1 10
8 1 1 1 -1 10
10 1 -1 -1 1 10
11 -1 1 -1 1 10
2 1 -1 -1 -1 11
3 -1 1 -1 -1 11
13 -1 -1 1 1 11
16 1 1 1 1 11
Use alias() to check for confounding with
blocks:
> alias(y ~ Block + A * B * C * D, a);
Model :
y ~ Block + A * B * C * D
Complete :
(Intercept) Block01 Block10 Block11 A B C D A:B A:C B:C A:D B:D A:C:D
C:D 1 -2 -2 0 0 0 0 0 0 0 0 0 0 0
A:B:C -1 0 2 2 0 0 0 0 0 0 0 0 0 0
A:B:D -1 2 0 2 0 0 0 0 0 0 0 0 0 0
B:C:D A:B:C:D
C:D 0 0
A:B:C 0 0
A:B:D 0 0
The numbers in the table give the details of the confounding; the important point is that all three confounded effects are listed, and are confounded only with the various block effects.
Partial Confounding
If blocks contain too few runs for the full basic design, but enough resources are available to replicate the blocked design, different effects may be confounded with blocks in each replicate; for example, a 23 design in 3 replicates, each in 2 blocks of 4 runs, with ABC confounded with blocks in Rep1, AB confounded with blocks in Rep2, and BC confounded with blocks in Rep3:
> rep1 = rep2 = rep3 = twoLevel(3);
> rep1$Block = paste(1, (1 + rep1$A * rep1$B * rep1$C) / 2, sep = "");
> rep2$Block = paste(2, (1 + rep1$A * rep1$B) / 2, sep = "");
> rep3$Block = paste(3, (1 + rep1$B * rep1$C) / 2, sep = "");
> a = rbind(rep1, rep2, rep3);
> a$Block = factor(a$Block);
> a[order(a$Block),]
A B C Block
1 -1 -1 -1 10
4 1 1 -1 10
6 1 -1 1 10
7 -1 1 1 10
2 1 -1 -1 11
3 -1 1 -1 11
5 -1 -1 1 11
8 1 1 1 11
10 1 -1 -1 20
11 -1 1 -1 20
14 1 -1 1 20
15 -1 1 1 20
9 -1 -1 -1 21
12 1 1 -1 21
13 -1 -1 1 21
16 1 1 1 21
19 -1 1 -1 30
20 1 1 -1 30
21 -1 -1 1 30
22 1 -1 1 30
17 -1 -1 -1 31
18 1 -1 -1 31
23 -1 1 1 31
24 1 1 1 31
Use alias() to check for confounding with
blocks:
> a$y = 0; > alias(y ~ Block + A * B * C, a); Model : y ~ Block + A * B * C
No effects are listed, but this means only that no effect is
completely confounded with blocks. To see the presence
of partial confounding, add the option partial =
TRUE:
> alias(y ~ Block + A * B * C, a, partial = TRUE);
Model :
y ~ Block + A * B * C
Partial :
(Intercept) Block11 Block20 Block21 Block30 Block31
(Intercept) 0 -0.7745967 -0.7071068 -0.7071068 -0.7071068 -0.7071068
Block11 0 0.0000000 0.5477226 0.5477226 0.5477226 0.5477226
Block20 0 0.0000000 0.0000000 0.4000000 0.5000000 0.5000000
Block21 0 0.0000000 0.0000000 0.0000000 0.5000000 0.5000000
Block30 0 0.0000000 0.0000000 0.0000000 0.0000000 0.4000000
Block31 0 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
A 0 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
B 0 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
C 0 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
A:B 0 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
A:C 0 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
B:C 0 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
A:B:C 0 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
A B C A:B A:C B:C A:B:C
(Intercept) 0 0 0 0.0000000 0 0.0000000 0.4472136
Block11 0 0 0 0.0000000 0 0.0000000 -0.5773503
Block20 0 0 0 0.3162278 0 0.0000000 -0.3162278
Block21 0 0 0 -0.3162278 0 0.0000000 -0.3162278
Block30 0 0 0 0.0000000 0 0.3162278 -0.3162278
Block31 0 0 0 0.0000000 0 -0.3162278 -0.3162278
A 0 0 0 0.0000000 0 0.0000000 0.0000000
B 0 0 0 0.0000000 0 0.0000000 0.0000000
C 0 0 0 0.0000000 0 0.0000000 0.0000000
A:B 0 0 0 0.0000000 0 0.0000000 0.0000000
A:C 0 0 0 0.0000000 0 0.0000000 0.0000000
B:C 0 0 0 0.0000000 0 0.0000000 0.0000000
A:B:C 0 0 0 0.0000000 0 0.0000000 0.0000000
The values in the upper right triangle of the table are correlations between estimated effects. The partially confounded effects (AB, BC, and ABC) are those that are correlated with the intercept or the block effects.
Fractional Factorial Designs
A 2k-p fractional factorial design is characterized by its full defining relation, but is usually constructed in steps:
- begin with a basic design in the first k-p factors;
- add columns for the remaining factors as interactions among the basic factors.
For example, a 29-5 Resolution-III design generated by E = BD, F = BCD, G = AC, H = ACD, J = AB:
> a = twoLevel(4); > a$E = a$B * a$D; > a$F = a$B * a$C * a$D; > a$G = a$A * a$C; > a$H = a$A * a$C * a$D; > a$J = a$A * a$B;
The alias() function shows the alias
structure:
> a$y = 0;
> alias(y ~ A * B * C * D * E * F * G * H * J, a)
Model :
y ~ A * B * C * D * E * F * G * H * J
Complete :
(Intercept) A B C D E F G H J B:C A:D C:D A:E A:F B:G
C:G 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0
D:G 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0
E:G 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0
F:G 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0
[many rows omitted]
A:B:D:E:F:G:H:J 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0
A:C:D:E:F:G:H:J 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0
B:C:D:E:F:G:H:J 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0
A:B:C:D:E:F:G:H:J 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
A:B 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0
A:C 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0
B:D 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0
B:E 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0
C:E 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0
D:E 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0
B:F 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0
C:F 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0
D:F 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0
E:F 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0
A:G 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0
The table has one column for each alias chain, labeled by
the effect that will be listed for instance in
summary(lm(y ~ A * B * C * D * E * F * G * H * J,
a)). It has one row for every other effect (496 such
rows in this example), which shows with which listed effect it
is aliased.
The alias chain for the intercept is the full defining relation of the design:
> aAlias = alias(y ~ A * B * C * D * E * F * G * H * J, a)$Complete; > aAlias = round(aAlias); > dimnames(aAlias)[[1]][aAlias[ ,"(Intercept)"] == 1]; [1] "B:D:E" "C:E:F" "A:C:G" [4] "D:G:H" "A:B:J" "F:H:J" [7] "B:C:D:F" "A:E:F:G" "A:C:D:H" [10] "A:B:F:H" "B:E:G:H" "A:D:E:J" [13] "B:C:G:J" "D:F:G:J" "C:E:H:J" [16] "A:B:D:F:G" "A:B:C:E:H" "A:D:E:F:H" [19] "B:C:F:G:H" "A:C:D:F:J" "C:D:E:G:J" [22] "B:E:F:G:J" "B:C:D:H:J" "A:E:G:H:J" [25] "A:B:C:D:E:G" "C:D:E:F:G:H" "A:B:C:E:F:J" [28] "B:D:E:F:H:J" "A:B:D:G:H:J" "A:C:F:G:H:J" [31] "A:B:C:D:E:F:G:H:J"
The alias chain for any effect can be extracted in the same way:
> dimnames(aAlias)[[1]][aAlias[, "A"] == 1]; [1] "C:G" "B:J" "E:F:G" "C:D:H" [5] "B:F:H" "D:E:J" "A:B:D:E" "A:C:E:F" [9] "B:D:F:G" "B:C:E:H" "D:E:F:H" "A:D:G:H" [13] "C:D:F:J" "A:F:H:J" "E:G:H:J" "A:B:C:D:F" [17] "B:C:D:E:G" "A:B:E:G:H" "B:C:E:F:J" "A:B:C:G:J" [21] "A:D:F:G:J" "A:C:E:H:J" "B:D:G:H:J" "C:F:G:H:J" [25] "A:B:C:F:G:H" "A:C:D:E:G:J" "A:B:E:F:G:J" "A:B:C:D:H:J" [29] "A:C:D:E:F:G:H" "A:B:D:E:F:H:J" "B:C:D:E:F:G:H:J"
In a Resolution-III design like this, aliasing is often presented in terms only of the two-factor interactions that are aliased with each main effect.
> aAliasA = dimnames(aAlias)[[1]][aAlias[, "A"] == 1];
> aAliasA[nchar(aAliasA) < 5];
[1] "C:G" "B:J"
