Assignment 16
IE 5342 | 2^k Factorial Designs using R
Jeffrey Mateo Hanna Hernandez
2026-11-04
Problem 1: Silver powder production
A silver powder process was run as a \(2^3\) factorial with two replicates. The response is powder density and the three factors are Ammonium (2% or 30%), Stir Rate (100 or 150 RPM) and Temperature (8 or 40 degrees C).
Data
dat <- read.csv("https://raw.githubusercontent.com/tmatis12/datafiles/main/PowderProduction.csv",
fileEncoding = "UTF-8-BOM")
A <- factor(dat$Ammonium)
B <- factor(dat$StirRate)
C <- factor(dat$Temperature)
dat <- data.frame(A, B, C, Density = dat$Density)
head(dat)## A B C Density
## 1 2 100 8 14.68
## 2 2 100 8 15.18
## 3 30 100 8 15.12
## 4 30 100 8 17.48
## 5 2 150 8 7.54
## 6 2 150 8 6.66str(dat)## 'data.frame': 16 obs. of 4 variables:
## $ A : Factor w/ 2 levels "2","30": 1 1 2 2 1 1 2 2 1 1 ...
## $ B : Factor w/ 2 levels "100","150": 1 1 1 1 2 2 2 2 1 1 ...
## $ C : Factor w/ 2 levels "8","40": 1 1 1 1 1 1 1 1 2 2 ...
## $ Density: num 14.68 15.18 15.12 17.48 7.54 ...The three columns are turned into factors because each one only has two levels and we want R to treat them as groups, not as numbers.
Mean density at each level:
tapply(dat$Density, dat$A, mean)## 2 30
## 11.19500 14.52625tapply(dat$Density, dat$B, mean)## 100 150
## 14.96250 10.75875tapply(dat$Density, dat$C, mean)## 8 40
## 12.71750 13.00375a. Model equation
The full \(2^3\) factorial model with two replicates is
\[ y_{ijkl} = \mu + \tau_i + \beta_j + \gamma_k + (\tau\beta)_{ij} + (\tau\gamma)_{ik} + (\beta\gamma)_{jk} + (\tau\beta\gamma)_{ijk} + \varepsilon_{ijkl} \]
where
- \(y_{ijkl}\) is the powder density of replicate \(l\)
- \(\mu\) is the overall mean
- \(\tau_i\) is the effect of Ammonium, \(i = 1, 2\)
- \(\beta_j\) is the effect of Stir Rate, \(j = 1, 2\)
- \(\gamma_k\) is the effect of Temperature, \(k = 1, 2\)
- \((\tau\beta)_{ij}\), \((\tau\gamma)_{ik}\) and \((\beta\gamma)_{jk}\) are the two factor interactions
- \((\tau\beta\gamma)_{ijk}\) is the three factor interaction
- \(\varepsilon_{ijkl}\) is the error, with \(\varepsilon_{ijkl} \sim N(0, \sigma^2)\) and \(l = 1, 2\)
All three factors are fixed.
b. Which factors are significant?
Full model
A*B*C tells R to fit the three main effects and all the
interactions at once.
full_model <- aov(Density ~ A * B * C, data = dat)
summary(full_model)## Df Sum Sq Mean Sq F value Pr(>F)
## A 1 44.39 44.39 11.180 0.01018 *
## B 1 70.69 70.69 17.804 0.00292 **
## C 1 0.33 0.33 0.083 0.78117
## A:B 1 28.12 28.12 7.082 0.02875 *
## A:C 1 0.02 0.02 0.005 0.94281
## B:C 1 10.13 10.13 2.551 0.14889
## A:B:C 1 1.52 1.52 0.383 0.55341
## Residuals 8 31.76 3.97
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Using alpha = 0.05, the terms that come out significant are:
| Term | p-value | Significant? |
|---|---|---|
| A (Ammonium) | 0.0102 | Yes |
| B (Stir Rate) | 0.0029 | Yes |
| C (Temperature) | 0.7812 | No |
| A:B | 0.0288 | Yes |
| A:C | 0.9428 | No |
| B:C | 0.1489 | No |
| A:B:C | 0.5534 | No |
Interaction plots
interaction.plot(dat$A, dat$B, dat$Density,
main = "Ammonium by Stir Rate",
xlab = "Ammonium (%)",
ylab = "Mean Density",
trace.label = "Stir Rate")
interaction.plot(dat$A, dat$C, dat$Density,
main = "Ammonium by Temperature",
xlab = "Ammonium (%)",
ylab = "Mean Density",
trace.label = "Temperature")
interaction.plot(dat$B, dat$C, dat$Density,
main = "Stir Rate by Temperature",
xlab = "Stir Rate (RPM)",
ylab = "Mean Density",
trace.label = "Temperature")
In the first plot the two lines are far from parallel. At 100 RPM the line is almost flat, so changing the ammonium barely moves the density. At 150 RPM the line goes up steeply, from 7.77 to 13.75. That is the A:B interaction, and it is the one the table marks as significant.
In the second plot the two lines are almost parallel, so there is no interaction between ammonium and temperature.
In the third plot the lines cross a little, but the p-value is 0.1489, so that small crossing is not enough to call it an interaction.
Which terms are dropped
We drop C, A:C, B:C and A:B:C because all of them have p-values above 0.05 (0.7812, 0.9428, 0.1489 and 0.5534). Temperature is the only main effect we drop. Its mean density is 12.72 at 8 degrees and 13.00 at 40 degrees, a difference of only 0.29, which is tiny next to the variation in the data.
We keep A, B and A:B. Note that even if a main effect came out weak, we would still have to keep A and B here, because A:B is significant and a model should not have an interaction without the main effects that form it.
Reduced model
reduced_model <- aov(Density ~ A * B, data = dat)
summary(reduced_model)## Df Sum Sq Mean Sq F value Pr(>F)
## A 1 44.39 44.39 12.17 0.004472 **
## B 1 70.69 70.69 19.38 0.000861 ***
## A:B 1 28.12 28.12 7.71 0.016751 *
## Residuals 12 43.76 3.65
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Model check
qqnorm(residuals(reduced_model), main = "Normal Q-Q Plot of the Residuals")
qqline(residuals(reduced_model))
plot(fitted(reduced_model), residuals(reduced_model),
main = "Residuals vs Fitted Values",
xlab = "Fitted Value",
ylab = "Residual",
pch = 19)
abline(h = 0, lty = 2)
The Q-Q plot is close to the line and the residuals vs fitted plot has no pattern and no funnel, so the reduced model is fine.
Final answer
The significant terms and their final p-values from the reduced model are:
| Term | Final p-value |
|---|---|
| A (Ammonium) | 0.004472 |
| B (Stir Rate) | 0.000861 |
| A:B | 0.016751 |
Density goes up when the ammonium goes from 2% to 30% and goes down when the stir rate goes from 100 to 150 RPM. Because the A:B interaction is significant, these two cannot be read separately: how much the ammonium helps depends on the stir rate. Temperature does not affect the density in the range that was tested.
tapply(dat$Density, list(dat$A, dat$B), mean)## 100 150
## 2 14.6225 7.7675
## 30 15.3025 13.7500The best cell is 30% ammonium at 100 RPM, with a mean density of 15.30.
Problem 2: Greenhouse tomato design
Here we only build the layout, we do not analyze any data. The experiment is a \(2^3\) factorial with three replications in a randomized complete block design.
The factors are:
| Factor | Description | Level 1 | Level 2 |
|---|---|---|---|
| A | Fertilizer Type | Organic compost | Chemical fertilizer |
| B | Irrigation Frequency | Once per day | Twice per day |
| C | Temperature Setting | 22 C | 28 C |
library(agricolae)The function design.ab() builds the layout. The argument trt = c(2, 2, 2) says three factors with two levels each, r = 3 is the three replications, design = “rcbd” makes the replications into blocks, and the seed makes the randomization repeatable.
design <- design.ab(trt = c(2, 2, 2),
r = 3,
serie = 2,
design = "rcbd",
seed = 20261104)
design$parameters## $design
## [1] "factorial"
##
## $trt
## [1] "1 1 1" "1 1 2" "1 2 1" "1 2 2" "2 1 1" "2 1 2" "2 2 1" "2 2 2"
##
## $r
## [1] 3
##
## $serie
## [1] 2
##
## $seed
## [1] 20261104
##
## $kinds
## [1] "Super-Duper"
##
## [[7]]
## [1] TRUE
##
## $applied
## [1] "rcbd"Layout
layout <- design$book
layout## plots block A B C
## 1 101 1 1 2 2
## 2 102 1 2 1 2
## 3 103 1 1 1 1
## 4 104 1 2 1 1
## 5 105 1 2 2 2
## 6 106 1 1 1 2
## 7 107 1 1 2 1
## 8 108 1 2 2 1
## 9 109 2 1 1 2
## 10 110 2 2 1 1
## 11 111 2 1 2 2
## 12 112 2 2 2 1
## 13 113 2 2 1 2
## 14 114 2 2 2 2
## 15 115 2 1 2 1
## 16 116 2 1 1 1
## 17 117 3 1 2 1
## 18 118 3 1 1 1
## 19 119 3 2 2 2
## 20 120 3 1 1 2
## 21 121 3 1 2 2
## 22 122 3 2 1 2
## 23 123 3 2 2 1
## 24 124 3 2 1 1The table has 24 rows, which is the 8 treatment combinations times
the 3 blocks. plots is the plot number, block
is the greenhouse bench, and A, B and C are the levels of each factor.
Inside each block the 8 combinations appear once, in random order, which
is what a randomized complete block design means.
table(layout$block)##
## 1 2 3
## 8 8 8Layout with the real level names
The same table is easier to read with the names from the problem instead of 1 and 2.
named <- layout
levels(named$A) <- c("Organic compost", "Chemical fertilizer")
levels(named$B) <- c("Once per day", "Twice per day")
levels(named$C) <- c("22 C", "28 C")
head(named, 8)## plots block A B C
## 1 101 1 1 2 2
## 2 102 1 2 1 2
## 3 103 1 1 1 1
## 4 104 1 2 1 1
## 5 105 1 2 2 2
## 6 106 1 1 1 2
## 7 107 1 1 2 1
## 8 108 1 2 2 1Complete code
# ---- Problem 1 ----
dat <- read.csv("https://raw.githubusercontent.com/tmatis12/datafiles/main/PowderProduction.csv",
fileEncoding = "UTF-8-BOM")
A <- factor(dat$Ammonium)
B <- factor(dat$StirRate)
C <- factor(dat$Temperature)
dat <- data.frame(A, B, C, Density = dat$Density)
head(dat)
str(dat)
tapply(dat$Density, dat$A, mean)
tapply(dat$Density, dat$B, mean)
tapply(dat$Density, dat$C, mean)
# full model with all interactions
full_model <- aov(Density ~ A * B * C, data = dat)
summary(full_model)
# interaction plots
interaction.plot(dat$A, dat$B, dat$Density,
main = "Ammonium by Stir Rate",
xlab = "Ammonium (%)",
ylab = "Mean Density",
trace.label = "Stir Rate")
interaction.plot(dat$A, dat$C, dat$Density,
main = "Ammonium by Temperature",
xlab = "Ammonium (%)",
ylab = "Mean Density",
trace.label = "Temperature")
interaction.plot(dat$B, dat$C, dat$Density,
main = "Stir Rate by Temperature",
xlab = "Stir Rate (RPM)",
ylab = "Mean Density",
trace.label = "Temperature")
# reduced model
reduced_model <- aov(Density ~ A * B, data = dat)
summary(reduced_model)
qqnorm(residuals(reduced_model), main = "Normal Q-Q Plot of the Residuals")
qqline(residuals(reduced_model))
plot(fitted(reduced_model), residuals(reduced_model),
main = "Residuals vs Fitted Values",
xlab = "Fitted Value",
ylab = "Residual",
pch = 19)
abline(h = 0, lty = 2)
tapply(dat$Density, list(dat$A, dat$B), mean)
# ---- Problem 2 ----
library(agricolae)
design <- design.ab(trt = c(2, 2, 2),
r = 3,
serie = 2,
design = "rcbd",
seed = 20261104)
design$parameters
layout <- design$book
layout
table(layout$block)
named <- layout
levels(named$A) <- c("Organic compost", "Chemical fertilizer")
levels(named$B) <- c("Once per day", "Twice per day")
levels(named$C) <- c("22 C", "28 C")
head(named, 8)