Assignment 14: Factorial Designs Using R
William Scroggins
Last compiled on October 28, 2026 at 8:53 PM - CDT
1 Problem 1
1.1 Data Import
The data will be imported using the read.csv() function in the code chunk below. Additionally, the factors will be converted to fixed effects.
dat1<-read.csv("https://raw.githubusercontent.com/tmatis12/datafiles/main/PowderProduction.csv")
#install.packages("GAD")
library(GAD)
dat1$Ammonium<-as.fixed(dat1$Ammonium)
dat1$StirRate<-as.fixed(dat1$StirRate)
dat1$Temperature<-as.fixed(dat1$Temperature)1.2 Model Equation
The model equation for this factorial design with three factors can be seen in the equation below.
\[ X_{ijkl}=\mu + \alpha_i + \beta_j + \gamma_k + (\alpha\beta)_{ij} + (\alpha\gamma)_{ik} + (\beta\gamma)_{jk} + (\alpha\beta\gamma)_{ijk} + \epsilon_{ijkl} \]
Where \(X_{ijkl}\) is the response variable density, \(\mu\) is the grand mean, \(\alpha_i\) is the effect of Ammonium content, \(\beta_j\) is the effect of the Stir Rate, \(\gamma_k\) is the effect of the Temperature, the different combinations of the main effects are the two and three factor interaction effects, and \(\epsilon_{ijkl}\) is the random error associated with the measurement.
1.3 Factorial Design Analysis
The analysis will be performed using the “GAD” package and can be seen in the code chunk below.
model1<-lm(Density~Ammonium + StirRate + Temperature + Ammonium:StirRate + Ammonium:Temperature + StirRate:Temperature + Ammonium:StirRate:Temperature, data = dat1)
gad(model1)## $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 ' ' 1The gad(model1) function returned a p-value of 0.553412 for the three factor interaction. Since the three factor interaction is not significant, it will be removed from the model before analyzing the two factor interactions. The updated model and analysis can be seen in the code chunk below.
model2<-lm(Density~Ammonium + StirRate + Temperature + Ammonium:StirRate + Ammonium:Temperature + StirRate:Temperature, data = dat1)
gad(model2)## $anova
## Analysis of Variance Table
##
## Response: Density
## Df Sum Sq Mean Sq F value Pr(>F)
## Ammonium 1 44.389 44.389 12.0037 0.007109 **
## StirRate 1 70.686 70.686 19.1150 0.001792 **
## Temperature 1 0.328 0.328 0.0886 0.772681
## Ammonium:StirRate 1 28.117 28.117 7.6033 0.022206 *
## Ammonium:Temperature 1 0.022 0.022 0.0059 0.940538
## StirRate:Temperature 1 10.128 10.128 2.7389 0.132317
## Residuals 9 33.281 3.698
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1The gad(model2) function returned a p-value of 0.022206 for the interaction effect between the Ammonium content and the Stir Rate. Since we are using a significance level of 0.05, this interaction between the two factors is deemed significant. An interaction plot will now be generated between the two interacting factors in the code chunk below.
interaction.plot(dat1$StirRate, dat1$Ammonium, dat1$Density, main = "Interaction Plot Between Ammonium Content and Stir Rate", xlab = "Stir Rate", ylab = "Density", trace.label = "Ammonium Content")
As seen in the plot, the slopes of the two Ammonium levels differ significantly. For 2% Ammonium, the Stir Rate significantly effects the measured density of the powder, whereas the Stir Rate for 30% Ammonium has a minimal effect comparatively. It appears a crossover point might occur around a Stir Rate of 90 RPM where the 2% Ammonium would actually have a higher density than the 30% Ammonium.
2 Problem 2
2.1 Data Input
The data given will need to be entered into a data frame to analyze it. The creation of this data frame can be seen in the code chunk below.
#install.packages("GAD")
library(GAD)
Position<-c(1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2)
Temp<-rep(seq(from = 800, to = 850, by = 25), 6)
Obs<-c(570, 1063, 565, 565, 1080, 510, 583, 1043, 590, 528, 988, 526, 547, 1026, 538, 521, 1004, 532)
dat2<-data.frame(Position, Temp, Obs)2.2 Fixed Effects Analysis
An analysis will be performed where we assume that the factors Temperature and Position are fixed effects. This analysis can be seen in the code chunk below.
dat2$Position<-as.fixed(dat2$Position)
dat2$Temp<-as.fixed(dat2$Temp)
model3<-lm(Obs~Position + Temp + Position:Temp, data = dat2)
gad(model3)## $anova
## Analysis of Variance Table
##
## Response: Obs
## Df Sum Sq Mean Sq F value Pr(>F)
## Position 1 7160 7160 15.998 0.001762 **
## Temp 2 945342 472671 1056.117 3.25e-14 ***
## Position:Temp 2 818 409 0.914 0.427110
## Residuals 12 5371 448
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1The gad(model3) function returned a p-value of 0.427110 for the interaction effect between the Furnace Position and the Baking Temperature. Since the two factor interaction is not significant, it will be removed from the model before analyzing the main effects. The updated model and analysis can be seen in the code chunk below.
model4<-lm(Obs~Position + Temp, data = dat2)
gad(model4)## $anova
## Analysis of Variance Table
##
## Response: Obs
## Df Sum Sq Mean Sq F value Pr(>F)
## Position 1 7160 7160 16.197 0.001254 **
## Temp 2 945342 472671 1069.257 4.924e-16 ***
## Residuals 14 6189 442
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1The gad(model4) function returned p-values of 0.001254 for the Position effect and 4.924e-16 for the Temperature effect. Both of these main effects are significant at a significance level of 0.05.
2.3 Random Effects Analysis
Now an analysis will be performed where we assume that the factors Temperature and Position are random effects. This analysis can be seen in the code chunk below.
dat2$Position<-as.random(dat2$Position)
dat2$Temp<-as.random(dat2$Temp)
model5<-lm(Obs~Position + Temp + Position:Temp, data = dat2)
gad(model5)## $anova
## Analysis of Variance Table
##
## Response: Obs
## Df Sum Sq Mean Sq F value Pr(>F)
## Position 1 7160 7160 17.504 0.0526583 .
## Temp 2 945342 472671 1155.518 0.0008647 ***
## Position:Temp 2 818 409 0.914 0.4271101
## Residuals 12 5371 448
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1The gad(model5) function returned a p-value of 0.4271101 for the interaction effect between the Furnace Position and the Baking Temperature. Since the two factor interaction is not significant, it will be removed from the model before analyzing the main effects. The updated model and analysis can be seen in the code chunk below.
model6<-lm(Obs~Position + Temp, data = dat2)
gad(model6)## $anova
## Analysis of Variance Table
##
## Response: Obs
## Df Sum Sq Mean Sq F value Pr(>F)
## Position 1 7160 7160 16.197 0.001254 **
## Temp 2 945342 472671 1069.257 4.924e-16 ***
## Residuals 14 6189 442
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1The gad(model6) function returned p-values of 0.001254 for the Position effect and 4.924e-16 for the Temperature effect. Both of these main effects are significant at a significance level of 0.05.
2.4 Mixed Effects Analysis
Lastly an analysis will be performed where we assume that Temperature is a random effect and Position is a fixed effect. This analysis can be seen in the code chunk below.
dat2$Position<-as.fixed(dat2$Position)
dat2$Temp<-as.random(dat2$Temp)
model7<-lm(Obs~Position + Temp + Position:Temp, data = dat2)
gad(model7)## $anova
## Analysis of Variance Table
##
## Response: Obs
## Df Sum Sq Mean Sq F value Pr(>F)
## Position 1 7160 7160 17.504 0.05266 .
## Temp 2 945342 472671 1056.117 3.25e-14 ***
## Position:Temp 2 818 409 0.914 0.42711
## Residuals 12 5371 448
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1The gad(model7) function returned a p-value of 0.42711 for the interaction effect between the Furnace Position and the Baking Temperature. Since the two factor interaction is not significant, it will be removed from the model before analyzing the main effects. The updated model and analysis can be seen in the code chunk below.
model8<-lm(Obs~Position + Temp, data = dat2)
gad(model8)## $anova
## Analysis of Variance Table
##
## Response: Obs
## Df Sum Sq Mean Sq F value Pr(>F)
## Position 1 7160 7160 16.197 0.001254 **
## Temp 2 945342 472671 1069.257 4.924e-16 ***
## Residuals 14 6189 442
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1The gad(model8) function returned p-values of 0.001254 for the Position effect and 4.924e-16 for the Temperature effect. Both of these main effects are significant at a significance level of 0.05.
2.5 Analyses Comparisons
The p-values for the two factor interaction effect stayed constant regardless of the factors being fixed or random. Additionally, the p-values for the main effects after the two factor interaction was removed from the model stayed constant regardless of the factors being fixed or random. The only instance in which the fixed or random effect changed the p-value was for the main effects when the two factor interaction was still included in the model.
3 Complete R Code
# Problem 1
## Data Import
dat1<-read.csv("https://raw.githubusercontent.com/tmatis12/datafiles/main/PowderProduction.csv")
#install.packages("GAD")
library(GAD)
dat1$Ammonium<-as.fixed(dat1$Ammonium)
dat1$StirRate<-as.fixed(dat1$StirRate)
dat1$Temperature<-as.fixed(dat1$Temperature)
## Factorial Design Analysis
model1<-lm(Density~Ammonium + StirRate + Temperature + Ammonium:StirRate + Ammonium:Temperature + StirRate:Temperature + Ammonium:StirRate:Temperature, data = dat1)
gad(model1)
model2<-lm(Density~Ammonium + StirRate + Temperature + Ammonium:StirRate + Ammonium:Temperature + StirRate:Temperature, data = dat1)
gad(model2)
interaction.plot(dat1$StirRate, dat1$Ammonium, dat1$Density, main = "Interaction Plot Between Ammonium Content and Stir Rate", xlab = "Stir Rate", ylab = "Density", trace.label = "Ammonium Content")
# Problem 2
## Data Input
#install.packages("GAD")
library(GAD)
Position<-c(1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2)
Temp<-rep(seq(from = 800, to = 850, by = 25), 6)
Obs<-c(570, 1063, 565, 565, 1080, 510, 583, 1043, 590, 528, 988, 526, 547, 1026, 538, 521, 1004, 532)
dat2<-data.frame(Position, Temp, Obs)
## Fixed Effects Analysis
dat2$Position<-as.fixed(dat2$Position)
dat2$Temp<-as.fixed(dat2$Temp)
model3<-lm(Obs~Position + Temp + Position:Temp, data = dat2)
gad(model3)
model4<-lm(Obs~Position + Temp, data = dat2)
gad(model4)
## Random Effects Analysis
dat2$Position<-as.random(dat2$Position)
dat2$Temp<-as.random(dat2$Temp)
model5<-lm(Obs~Position + Temp + Position:Temp, data = dat2)
gad(model5)
model6<-lm(Obs~Position + Temp, data = dat2)
gad(model6)
## Mixed Effects Analysis
dat2$Position<-as.fixed(dat2$Position)
dat2$Temp<-as.random(dat2$Temp)
model7<-lm(Obs~Position + Temp + Position:Temp, data = dat2)
gad(model7)
model8<-lm(Obs~Position + Temp, data = dat2)
gad(model8)