Factorial Experiments
"Factorial experiment" is a general term used for experiments with multiple factors, designed so that all possible combinations of levels of factors are tested. In theory, such experiments could be quite large. The number of possible combinations of levels is given by the product of the number of levels in each factor. For example if
then the total number of combinations of levels is .
Note that this number grows quite quickly... Even if only
Observations are called replicates if they involve the same combination of levels/treatments of each factor.
Notation
- Upper case letters for factors,
e.g. , , when there are few factors
and when there are many factors - Lower case letters for the number of levels/treatments in each factor
e.g. or - Numbers for individual levels/treatments in factors
e.g. , or etc.
If there are only two levels, then we may use different notation like or else (see the section on factorial experiments) total number of observations number of replicate observations for each combination of levels.
So . will be indices for factors (unless we use for total number of factors)- Observations are written as follows (for a three factor experiment):
replicate observation with factor levels , , and
i.e. first factor has level , second factor has level , etc.
- Sample means are as follows (for experiment with factors
, , ): sample mean of all replicates with , ,
This is the mean of the cell sample mean of all observations with , and
This is the mean of the line sample mean of all observations with
This is the mean of the plane- etc.
We mentally organize data into a
Interaction between factors
If the effects of changing two factors are not independent, then we say there is interaction between them. Factorial experiments are designed to look for interaction between factors.
The factorial model assumes that expected values for each combination of factor values can be written as the overall mean, plus a main (treatment) effect of changing individual factors, plus interaction effects of changing combinations of factors.
Our general model looks something like the following:
- "Model:"
is the over-all mean are the "main treatment effects" for the different levels of corresponding factors are the "interaction effects" of the corresponding factors is the error, distributed as .
If two factors have no interaction, then when we plot their sample means we should get a series of parallel shapes. If there is interaction, then slopes will be unequal.
### Example: two factors
A <- c(1,2,3) # Levels of factor A
B <- c(1,2,3,4) # Levels of factor B
a <- length(A)
b <- length(B)
## choose mean and treatment effects
mu <- runif(1, -2,2) # mean
tA <- runif(a, -4,4) # treatment A effects
tB <- runif(b, -4,4) # treatment B effects
## combine mean and treatment effects for each level
Ybar <- mu +
matrix(rep(tA,b), nrow=a, byrow=FALSE) +
matrix(rep(tB,a), nrow=a, byrow=TRUE)
cat('means of experiment (no interaction)')
Ybar
## make factor indicators for each data point
fA <- matrix(rep(A,b), nrow=a, byrow=FALSE)
fB <- matrix(rep(B,a), nrow=a, byrow=TRUE)
## check that indicators are correct :)
cat('factor A indicator')
fA
cat('factor B indicator')
fB
## plot resulting means (no interaction)
plot(fA,Ybar,main="No interaction (parallel)")
## use lines to show B factor value
for (j in 1:b) {
lines(fA[,j],Ybar[,j], lty=j)
}
## now let's add an interaction term!
tAB <- matrix(runif(a*b, -4,4), nrow=a)
Ybar <- Ybar + tAB
## plot resulting means (with interaction)
plot(fA,Ybar,main="With interaction (not parallel)")
## use lines to show B factor value
for (j in 1:b) {
lines(fA[,j],Ybar[,j], lty=j)
}
Two Factor Experiments
Example: Observations for 2 factor experiment
Factor | B | |||
---|---|---|---|---|
Level | 1 | 2 | 3 | |
A | 1 | (1,1) observations | (1,2) observations | (1,3) observations |
2 | (2,1) observations | (2,2) observations | (2,3) observations | |
3 | (3,1) observations | (3,2) observations | (3,3) observations |
where "
For example, the upper-right corner of observations for
- The mean of each box is its replicate mean
- The means of the rows are the means of
observations: , , - The means of the columns are the means of
observations: , , - The overall mean is
Means | Means | |||
---|---|---|---|---|
Means |
(This table should remind you of the marginal probability tables that we made for joint random variables!)
Model
where
Alternately, the
where
Point estimators
-- the overall mean. -- the main treatment effect at -- the main treatment effect at -- the interaction effect at and .
Comparison to "single factor with blocking"
Note that a two factor experiment is very similar to a single factor experiment with blocking. There are two main differences, however.
- single factor with blocking experiments had no replicates -- within each level, each observation had a unique block.
- single factor with blocking experiments were assumed to have "level" and "block" effects independent -- there was no interaction effect
Sum of Squares
Setting variances equal on either side of the model equation leads to the "sum of square" terms.
- From
we get variance of Factor A (row) means:
- From
we get variance of Factor B (column) means: - From
we get variance of box means not accounted for by Factor A and B (row and column) variances: - From
we get pooled variance within boxes: - On the other side of the equation,
gives the variance of all observations:
Sources of Variation
Source | df | Sum of Square | Mean Square |
---|---|---|---|
Factor A | |||
Factor B | |||
Interaction | |||
Error | |||
Total |
(*) Not used.
Note:
F Tests
There are three different
- Model:
There are two tests for the treatment effects.
-
H0:
for all -
H0:
for all
And one test for interaction effect.
- H0:
for all
### Example: two factors (random experiment)
A <- c(1,2,3) # Levels of factor A
B <- c(1,2,3,4) # Levels of factor B
a <- length(A)
b <- length(B)
## choose random mean and treatment effects
mu <- runif(1, -2,2) # overall mean
tA <- runif(a, -2,2) # treatment A effects
tB <- runif(b, -2,2) # treatment B effects
tAB <- runif(a*b, -2,2) # interaction effects
## combine mean and treatment effects for each level
Ybar <- mu +
matrix(rep(tA,b), nrow=a, byrow=FALSE) +
matrix(rep(tB,a), nrow=a, byrow=TRUE) +
matrix(tAB, nrow=a)
## make factor indicators for each data point
fA <- matrix(rep(A,b), nrow=a, byrow=FALSE)
fB <- matrix(rep(B,a), nrow=a, byrow=TRUE)
## now add replicates for each level combination
n <- 50
y <- as.vector(rep(Ybar, n)) + rnorm(a*b*n, 0,5) # add error
fA <- factor( rep(fA, n))
fB <- factor( rep(fB, n))
data.AB <- data.frame(obs=y, A=fA, B=fB)
head(data.AB,10)
Example: Testing for interactions by hand!
### Let's compute the source of variance table by hand!!!!
ybar_A <- aggregate(obs ~ A, data=data.AB, mean)
ybar_B <- aggregate(obs ~ B, data=data.AB, mean)
ybar_AB <- aggregate(obs ~ A + B, data=data.AB, mean)
## sum of squares A and B are just var of row and col means
SS_A <- n*b * (a-1) * sd( ybar_A$obs )^2
SS_B <- n*a * (b-1) * sd( ybar_B$obs )^2
## sum of square interaction is "unexplained" variance of box means
SS_AB <- n * (a*b-1)* sd( ybar_AB$obs )^2 - SS_A - SS_B
## total sum of squares is just variance of data
SS_T <- (a*b*n-1) * sd( y )^2
## sum of square error is the "unexplained" part of total
SS_E <- SS_T - SS_A - SS_B - SS_AB
## convert to mean square
MS_A <- SS_A / (a-1)
MS_B <- SS_B / (b-1)
MS_AB <- SS_AB / ((a-1)*(b-1))
MS_E <- SS_E / (a*b*(n-1))
MS_T <- SS_T / (a*b*n-1)
## F-statistics and p-values!
F_A <- MS_A / MS_E
F_B <- MS_B / MS_E
F_AB <- MS_AB / MS_E
p_A <- 1 - pf(F_A, a-1, a*b*(n-1))
p_B <- 1 - pf(F_B, b-1, a*b*(n-1))
p_AB <- 1 - pf(F_AB, (a-1)*(b-1), a*b*(n-1))
## pretty table!
data.frame(
Source=c("Factor A","Factor B","Interaction", "Error", "Total"),
df =c( a-1, b-1, (a-1)*(b-1), a*b*(n-1),a*b*n-1),
SS =c( SS_A, SS_B, SS_AB, SS_E, SS_T ),
MS =c( MS_A, MS_B, MS_AB, MS_E, MS_T ),
F =c( F_A, F_B, F_AB, '', '' ),
p =c( p_A, p_B, p_AB, '', '' )
)
Example: Testing for interactions using R commands
## This is very quick and easy to do in R...
## interactions terms are indicated by multiplication
summary( aov(obs ~ A + B + A*B, data=data.AB) )
Example: Testing for interactions that don't exist
If interactions don't exist then we should be able to see this from the corresponding
An interesting thing happens when we remove treatment or interaction effects from our model. The corresponding row of the "Sources of Variance" table is removed and added to the Error row, without changing any of the other rows! This is called model collapse.
## What if there are no interactions??
### Example: two factors (random experiment)
A <- c(1,2,3) # Levels of factor A
B <- c(1,2,3,4) # Levels of factor B
a <- length(A)
b <- length(B)
## choose random mean and treatment effects
mu <- runif(1, -2,2) # overall mean
tA <- runif(a, -2,2) # treatment A effects
tB <- runif(b, -2,2) # treatment B effects
tAB <- rep(0, a*b ) # NO INTERACTION EFFECTS!!!!
## combine mean and treatment effects for each level
Ybar <- mu +
matrix(rep(tA,b), nrow=a, byrow=FALSE) +
matrix(rep(tB,a), nrow=a, byrow=TRUE) +
matrix(tAB, nrow=a)
## make factor indicators for each data point
fA <- matrix(rep(A,b), nrow=a, byrow=FALSE)
fB <- matrix(rep(B,a), nrow=a, byrow=TRUE)
## now add replicates for each level combination
n <- 50
y <- as.vector(rep(Ybar, n)) + rnorm(a*b*n, 0,5) # add error
fA <- factor( rep(fA, n))
fB <- factor( rep(fB, n))
data.noAB <- data.frame(obs=y, A=fA, B=fB)
cat('aov test with interactions')
summary( aov(obs ~ A + B + A*B, data=data.noAB) )
cat('aov test without interactions')
summary( aov(obs ~ A + B, data=data.noAB) )
Two factor experiments with no replicates ("one observation per cell")
- This is also called an "unreplicated experiment"
DO NOT DO THIS with only two factors!!!!!
In this case, the degrees of freedom of the error variance
If you assume that there is no interaction effect, then the interaction row will collapse onto the error row, so that we get some degrees of freedom. This is equivalent to one factor with blocking:
Model
Sources of Variation (no replicates, no interaction)
Source | df | Sum of Square | Mean Square |
---|---|---|---|
Factor A | |||
Factor B | |||
Error | |||
Total |
(*) if
... but one of the main points of factorial experiments is to look for interaction effects. So this is stupid.
Summary:
Without replicates, you cannot check for interaction effects in two factor experiments!
Three Factor Experiments
... in this case, data is a 3d array of lists of replicates.
Model
As well as the three main treatment effects, there are now three pair-wise interaction effects as well as one triple interaction effect.
Note that if there are
Point estimators
Point estimators of main effects and interactions are given by alternating sums:
Sum of Squares
Setting variances equal on either side of the model equation leads to the "sum of square" terms from the point estimators. For example:
Sources of Variation
Source | df | Sum of Square | Mean Square |
---|---|---|---|
Factor A | |||
Factor B | |||
Factor C | |||
Mix AB | |||
Mix BC | |||
Mix AC | |||
Mix ABC | |||
Error | |||
Total |
(*) Not used.
Now there are seven
Example. Random three factor experiment
First I'll generate treatment effects and data.
### Example: three factors (random experiment)
A <- c(1,2,3) # Levels of factor A
B <- c(1,2,3,4) # Levels of factor B
C <- c(1,2,3,4) # Levels of factor C
a <- length(A)
b <- length(B)
c <- length(C)
## choose random mean and treatment effects
mu <- runif(1, -2,2) # overall mean
tA <- runif(a, -2,2) # treatment A main effects
tB <- runif(b, -2,2) # treatment B main effects
tC <- runif(c, -2,2) # treatment C main effects
tAB <- runif(a*b, -2,2) # AB interaction effects
tAC <- runif(a*c, -2,2) # AC interaction effects
tBC <- runif(b*c, -2,2) # BC interaction effects
tABC<- runif(a*b*c, -2,2) # ABC interaction effects
## combine mean and treatment effects
## use array(.., dim=c(..)) to arrange into 3d arrays
## use aperm(.., c(..)) to reorder dimensions of arrays
Ybar <- mu +
array(tA, dim=c(a,b,c)) +
aperm(array(tB, dim=c(b,c,a)), c(3,1,2)) +
aperm(array(tC, dim=c(c,a,b)), c(2,3,1)) +
array(tAB, dim=c(a,b,c)) +
aperm(array(tAC, dim=c(a,c,b)), c(1,3,2)) +
aperm(array(tBC, dim=c(b,c,a)), c(3,1,2)) +
array(tABC, dim=c(a,b,c))
## compare to simpler versions missing some interactions
Ybar.noABC <- mu +
array(tA, dim=c(a,b,c)) +
aperm(array(tB, dim=c(b,c,a)), c(3,1,2)) +
aperm(array(tC, dim=c(c,a,b)), c(2,3,1)) +
array(tAB, dim=c(a,b,c)) +
aperm(array(tAC, dim=c(a,c,b)), c(1,3,2)) +
aperm(array(tBC, dim=c(b,c,a)), c(3,1,2))
Ybar.noBC <- mu +
array(tA, dim=c(a,b,c)) +
aperm(array(tB, dim=c(b,c,a)), c(3,1,2)) +
aperm(array(tC, dim=c(c,a,b)), c(2,3,1)) +
array(tAB, dim=c(a,b,c)) +
aperm(array(tAC, dim=c(a,c,b)), c(1,3,2))
Ybar.noB <- mu +
array(tA, dim=c(a,b,c)) +
aperm(array(tC, dim=c(c,a,b)), c(2,3,1)) +
aperm(array(tAC, dim=c(a,c,b)), c(1,3,2))
Ybar.onlyAC <- mu +
array(tA, dim=c(a,b,c)) +
aperm(array(tB, dim=c(b,c,a)), c(3,1,2)) +
aperm(array(tC, dim=c(c,a,b)), c(2,3,1)) +
aperm(array(tAC, dim=c(a,c,b)), c(1,3,2))
## make factor indicators for each data point
fA <- array(A, dim=c(a,b,c))
fB <- aperm(array(B, dim=c(b,c,a)), c(3,1,2))
fC <- aperm(array(C, dim=c(c,a,b)), c(2,3,1))
## now add replicates
n <- 50
## random error for replicates
error <- rnorm(a*b*c*n, 0,5)
## observations
y <- as.vector(rep(Ybar, n)) + error
## also add replicates for reduced versions
y.noABC <- as.vector(rep(Ybar.noABC, n)) + error
y.noBC <- as.vector(rep(Ybar.noBC, n)) + error
y.noB <- as.vector(rep(Ybar.noB, n)) + error
y.onlyAC <- as.vector(rep(Ybar.onlyAC, n)) + error
## now add replicates for factors
fA <- factor(rep(fA, n))
fB <- factor(rep(fB, n))
fC <- factor(rep(fC, n))
## combine everything into a big data frame!
data.ABC <- data.frame(
obs = y,
noABC = y.noABC,
noBC = y.noBC,
noB = y.noB,
onlyAC = y.onlyAC,
A=fA, B=fB, C=fC
)
Now let's look at some anova tables!