Assignment 14 - Factorial Designs
G M Farhan Tanvir
2026-10-28
Problem 1: Powder Production
(a) Full Factorial Model
The experiment studies the effects of ammonium concentration, stir rate, and temperature on powder density. Each factor has two levels, all factors are fixed effects, and each treatment combination is replicated twice.
The full factorial model is
\[ \begin{aligned} Y_{ijkl}={}&\mu+\alpha_i+\beta_j+\gamma_k\\ &+(\alpha\beta)_{ij}+(\alpha\gamma)_{ik}\\ &+(\beta\gamma)_{jk}+(\alpha\beta\gamma)_{ijk}\\ &+\epsilon_{ijkl}. \end{aligned} \]
Here, \(\mu\) is the overall mean; \(\alpha_i,\beta_j,\gamma_k\) are the main effects; the remaining terms represent interactions; and \(\epsilon_{ijkl}\) is random error.
Enter the Data
powder <- data.frame(
Ammonium = factor(c(
2, 2, 30, 30, 2, 2, 30, 30,
2, 2, 30, 30, 2, 2, 30, 30
)),
StirRate = factor(c(
100, 100, 100, 100, 150, 150, 150, 150,
100, 100, 100, 100, 150, 150, 150, 150
)),
Temperature = factor(c(
rep(8, 8), rep(40, 8)
)),
Density = c(
14.68, 15.18, 15.12, 17.48,
7.54, 6.66, 12.46, 12.62,
10.95, 17.68, 12.65, 15.96,
8.03, 8.84, 14.96, 14.96
)
)
powder## Ammonium StirRate Temperature 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.66
## 7 30 150 8 12.46
## 8 30 150 8 12.62
## 9 2 100 40 10.95
## 10 2 100 40 17.68
## 11 30 100 40 12.65
## 12 30 100 40 15.96
## 13 2 150 40 8.03
## 14 2 150 40 8.84
## 15 30 150 40 14.96
## 16 30 150 40 14.96(b) Factorial ANOVA
The significance level is \(\alpha=0.05\).
For each main effect and interaction:
\[ H_0:\text{The corresponding effect is zero} \]
\[ H_a:\text{The corresponding effect is not zero} \]
model1 <- aov(
Density ~ Ammonium * StirRate * Temperature,
data = powder
)
anova1 <- anova(model1)
anova1## Analysis of Variance Table
##
## Response: Density
## Df Sum Sq Mean Sq F value Pr(>F)
## Ammonium 1 44.389 44.389 11.1803 0.010175 *
## StirRate 1 70.686 70.686 17.8037 0.002918 **
## Temperature 1 0.328 0.328 0.0826 0.781170
## Ammonium:StirRate 1 28.117 28.117 7.0817 0.028754 *
## Ammonium:Temperature 1 0.022 0.022 0.0055 0.942808
## StirRate:Temperature 1 10.128 10.128 2.5510 0.148890
## Ammonium:StirRate:Temperature 1 1.519 1.519 0.3826 0.553412
## Residuals 8 31.762 3.970
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Significant Effects
results1 <- data.frame(
Effect = rownames(anova1),
F_value = anova1[["F value"]],
P_value = anova1[["Pr(>F)"]]
)
print(results1, row.names = FALSE)## Effect F_value P_value
## Ammonium 11.180257443 0.010175227
## StirRate 17.803734610 0.002918021
## Temperature 0.082552141 0.781170009
## Ammonium:StirRate 7.081719394 0.028754147
## Ammonium:Temperature 0.005479758 0.942807665
## StirRate:Temperature 2.551021886 0.148889937
## Ammonium:StirRate:Temperature 0.382605506 0.553411659
## Residuals NA NAsignificant1 <- subset(
results1,
!is.na(P_value) & P_value < 0.05
)
cat("\nSignificant effects at alpha = 0.05:\n")##
## Significant effects at alpha = 0.05:if (nrow(significant1) > 0) {
print(significant1, row.names = FALSE)
} else {
cat("No effects are significant at alpha = 0.05.\n")
}## Effect F_value P_value
## Ammonium 11.180257 0.010175227
## StirRate 17.803735 0.002918021
## Ammonium:StirRate 7.081719 0.028754147An effect is statistically significant when its p-value is below 0.05. The table identifies significant main effects and interactions affecting powder density.
Interaction Plots
par(mfrow = c(1, 3))
interaction.plot(
x.factor = powder$StirRate,
trace.factor = powder$Ammonium,
response = powder$Density,
fun = mean,
xlab = "Stir Rate (RPM)",
ylab = "Mean Density",
trace.label = "Ammonium (%)",
main = "Ammonium x Stir Rate"
)
interaction.plot(
x.factor = powder$Temperature,
trace.factor = powder$Ammonium,
response = powder$Density,
fun = mean,
xlab = "Temperature (C)",
ylab = "Mean Density",
trace.label = "Ammonium (%)",
main = "Ammonium x Temperature"
)
interaction.plot(
x.factor = powder$Temperature,
trace.factor = powder$StirRate,
response = powder$Density,
fun = mean,
xlab = "Temperature (C)",
ylab = "Mean Density",
trace.label = "Stir Rate",
main = "Stir Rate x Temperature"
)
par(mfrow = c(1, 1))Approximately parallel lines suggest little interaction. Nonparallel lines suggest a possible interaction. Statistical significance should be determined from the ANOVA p-values.
Problem 2: Carbon Anode Density
The experiment investigates whether firing temperature and furnace position affect baked carbon-anode density. There are two positions, three temperatures, and three observations per combination.
Prepare the Data
anode <- data.frame(
Position = factor(rep(c(1, 2), each = 9)),
Temperature = factor(
rep(rep(c(800, 825, 850), each = 3), times = 2)
),
Density = c(
# Position 1
570, 565, 583,
1063, 1080, 1043,
565, 510, 590,
# Position 2
528, 547, 521,
988, 1026, 1004,
526, 538, 532
)
)
anode## Position Temperature Density
## 1 1 800 570
## 2 1 800 565
## 3 1 800 583
## 4 1 825 1063
## 5 1 825 1080
## 6 1 825 1043
## 7 1 850 565
## 8 1 850 510
## 9 1 850 590
## 10 2 800 528
## 11 2 800 547
## 12 2 800 521
## 13 2 825 988
## 14 2 825 1026
## 15 2 825 1004
## 16 2 850 526
## 17 2 850 538
## 18 2 850 532The factorial model includes position, temperature, and their interaction:
\[ Y_{ijk}=\mu+\alpha_i+\beta_j+ (\alpha\beta)_{ij}+\epsilon_{ijk}. \]
(a) Both Factors Fixed
anode_a <- anode
anode_a$Position <- as.fixed(anode_a$Position)
anode_a$Temperature <- as.fixed(anode_a$Temperature)
model_a <- lm(
Density ~ Position * Temperature,
data = anode_a
)
result_a <- gad(model_a)
cat("Part (a): Both factors are fixed effects\n")## Part (a): Both factors are fixed effectsprint(result_a$anova)## Analysis of Variance Table
##
## Response: Density
## Df Sum Sq Mean Sq F value Pr(>F)
## Position 1 7160 7160 15.998 0.001762 **
## Temperature 2 945342 472671 1056.117 3.25e-14 ***
## Position:Temperature 2 818 409 0.914 0.427110
## Residuals 12 5371 448
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1(b) Both Factors Random
anode_b <- anode
anode_b$Position <- as.random(anode_b$Position)
anode_b$Temperature <- as.random(anode_b$Temperature)
model_b <- lm(
Density ~ Position * Temperature,
data = anode_b
)
result_b <- gad(model_b)
cat("Part (b): Both factors are random effects\n")## Part (b): Both factors are random effectsprint(result_b$anova)## Analysis of Variance Table
##
## Response: Density
## Df Sum Sq Mean Sq F value Pr(>F)
## Position 1 7160 7160 17.504 0.0526583 .
## Temperature 2 945342 472671 1155.518 0.0008647 ***
## Position:Temperature 2 818 409 0.914 0.4271101
## Residuals 12 5371 448
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1(c) Position Fixed and Temperature Random
anode_c <- anode
anode_c$Position <- as.fixed(anode_c$Position)
anode_c$Temperature <- as.random(anode_c$Temperature)
model_c <- lm(
Density ~ Position * Temperature,
data = anode_c
)
result_c <- gad(model_c)
cat("Part (c): Position fixed, Temperature random\n")## Part (c): Position fixed, Temperature randomprint(result_c$anova)## Analysis of Variance Table
##
## Response: Density
## Df Sum Sq Mean Sq F value Pr(>F)
## Position 1 7160 7160 17.504 0.05266 .
## Temperature 2 945342 472671 1056.117 3.25e-14 ***
## Position:Temperature 2 818 409 0.914 0.42711
## Residuals 12 5371 448
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1(d) Compare the P-values
The GAD output stores its ANOVA table in
$anova. The following code extracts the p-values directly
from that component.
tab_a <- result_a$anova
tab_b <- result_b$anova
tab_c <- result_c$anova
# Verify that all tables contain the expected p-value column.
p_column <- "Pr(>F)"
if (!p_column %in% names(tab_a) ||
!p_column %in% names(tab_b) ||
!p_column %in% names(tab_c)) {
stop("The expected p-value column was not found. Check the three GAD tables above.")
}
p_a <- setNames(tab_a[[p_column]], rownames(tab_a))
p_b <- setNames(tab_b[[p_column]], rownames(tab_b))
p_c <- setNames(tab_c[[p_column]], rownames(tab_c))
# Compare effects common to all three tables.
effects <- Reduce(
intersect,
list(names(p_a), names(p_b), names(p_c))
)
effects <- effects[
!grepl("Residuals", effects, ignore.case = TRUE)
]
comparison <- data.frame(
Effect = effects,
Both_Fixed = unname(p_a[effects]),
Both_Random = unname(p_b[effects]),
Position_Fixed_Temperature_Random = unname(p_c[effects])
)
cat("Comparison of p-values across the three models:\n")## Comparison of p-values across the three models:print(comparison, row.names = FALSE)## Effect Both_Fixed Both_Random
## Position 1.762434e-03 0.0526583234
## Temperature 3.249924e-14 0.0008646645
## Position:Temperature 4.271101e-01 0.4271101174
## Position_Fixed_Temperature_Random
## 5.265832e-02
## 3.249924e-14
## 4.271101e-01Discussion
The p-values are compared across the three models to assess how treating factors as fixed or random changes the conclusions.
Fixed effects concern the specific factor levels included in the experiment. Random effects concern variability among levels considered representative of a larger population. The mixed model treats position as fixed and temperature as random.
The comparison table reports p-values for the effects present in all three analyses. Interpret each p-value using the corresponding model and its assumptions.