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
\[y_{ijkl} = \mu + \alpha_i + \beta_j + \gamma_k + (\alpha\beta)_{ij} + (\alpha\gamma)_{ik} + (\beta\gamma)_{jk} + (\alpha\beta\gamma)_{ijk} + \varepsilon_{ijkl}, \qquad i, j, k, l = 1, 2\]
where \(\alpha_i\) = Ammonium, \(\beta_j\) = Stir Rate, \(\gamma_k\) = Temperature (all fixed), and \(\varepsilon_{ijkl} \sim NID(0, \sigma^2)\).
Full model
$anova
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 ' ' 1
The non-significant terms were removed one at a time (ABC → AC → BC → Temperature).
Reduced model
$anova
Analysis of Variance Table
Response: Density
Df Sum Sq Mean Sq F value Pr(>F)
Ammonium 1 44.389 44.389 12.1727 0.0044721 **
StirRate 1 70.686 70.686 19.3841 0.0008612 ***
Ammonium:StirRate 1 28.117 28.117 7.7103 0.0167511 *
Residuals 12 43.759 3.647
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
At α = 0.05, the significant factors are:
Temperature is not significant.
Increasing stir rate lowers density, but the drop is much larger at 2% ammonium than at 30% ammonium.
Temperature Rep Position Density
1 800 1 1 570
2 825 1 1 1063
3 850 1 1 565
4 800 2 1 565
5 825 2 1 1080
6 850 2 1 510
$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
Position: p = 0.0018, Temperature: p = 3.25 × 10⁻¹⁴, Interaction: p = 0.4271
$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
Position: p = 0.0527, Temperature: p = 0.00086, Interaction: p = 0.4271
$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
Position: p = 0.0527, Temperature: p = 3.25 × 10⁻¹⁴, Interaction: p = 0.4271
The interaction p-value (0.4271) is the same in all three cases because it is always tested against MS_E. The main effects change because the denominator of the F test changes:
Temperature is highly significant in every case, but the conclusion for Position depends on whether the effects are fixed or random.
library(GAD)
library(dplyr)
library(ggplot2)
# Question 1
dat <- read.csv("https://raw.githubusercontent.com/tmatis12/datafiles/main/PowderProduction.csv",
fileEncoding = "UTF-8-BOM")
dat$Ammonium <- as.fixed(dat$Ammonium)
dat$StirRate <- as.fixed(dat$StirRate)
dat$Temperature <- as.fixed(dat$Temperature)
head(dat)
# Full model
model_full <- lm(Density ~ Ammonium * StirRate * Temperature, data = dat)
gad(model_full)
# Reduced model
model_reduced <- lm(Density ~ Ammonium * StirRate, data = dat)
gad(model_reduced)
# Interaction plot
dat %>%
group_by(Ammonium, StirRate) %>%
summarise(Mean = mean(Density), .groups = "drop") %>%
ggplot(aes(x = StirRate, y = Mean, group = Ammonium, color = Ammonium)) +
geom_point(size = 2) +
geom_line() +
labs(x = "Stir Rate (RPM)", y = "Mean Density", color = "Ammonium (%)",
title = "Interaction Plot: Ammonium x Stir Rate")
# Question 2
dat2 <- expand.grid(Temperature = c(800, 825, 850), Rep = 1:3, Position = 1:2)
dat2$Density <- c(570, 1063, 565,
565, 1080, 510,
583, 1043, 590,
528, 988, 526,
547, 1026, 538,
521, 1004, 532)
head(dat2)
# a) Both fixed
dat_a <- dat2
dat_a$Position <- as.fixed(dat_a$Position)
dat_a$Temperature <- as.fixed(dat_a$Temperature)
model_a <- lm(Density ~ Position * Temperature, data = dat_a)
gad(model_a)
# b) Both random
dat_b <- dat2
dat_b$Position <- as.random(dat_b$Position)
dat_b$Temperature <- as.random(dat_b$Temperature)
model_b <- lm(Density ~ Position * Temperature, data = dat_b)
gad(model_b)
# c) Position fixed, Temperature random
dat_c <- dat2
dat_c$Position <- as.fixed(dat_c$Position)
dat_c$Temperature <- as.random(dat_c$Temperature)
model_c <- lm(Density ~ Position * Temperature, data = dat_c)
gad(model_c)