Factorial Designs using R

1.a. Let \(A\), \(B\), and \(C\) represent Ammonium, Stir Rate, and Temperature, respectively. The full factorial model is

\[ Y_{ijkl} = \mu +\alpha_i +\beta_j +\gamma_k +(\alpha\beta)_{ij} +(\alpha\gamma)_{ik} +(\beta\gamma)_{jk} +(\alpha\beta\gamma)_{ijk} +\epsilon_{ijkl}, \]

where \(i,j,k=1,2\), \(l=1,2\), and \(\epsilon_{ijkl}\sim N(0,\sigma^2)\).

1.b.

powder <- read.csv("https://raw.githubusercontent.com/tmatis12/datafiles/main/PowderProduction.csv")

powder$Ammonium <- factor(powder$Ammonium)
powder$StirRate <- factor(powder$StirRate)
powder$Temperature <- factor(powder$Temperature)

model <- aov(Density ~ Ammonium * StirRate * Temperature, data = powder)
summary(model)

##                               Df Sum Sq Mean Sq F value  Pr(>F)   
## Ammonium                       1  44.39   44.39  11.180 0.01018 * 
## StirRate                       1  70.69   70.69  17.804 0.00292 **
## Temperature                    1   0.33    0.33   0.083 0.78117   
## Ammonium:StirRate              1  28.12   28.12   7.082 0.02875 * 
## Ammonium:Temperature           1   0.02    0.02   0.005 0.94281   
## StirRate:Temperature           1  10.13   10.13   2.551 0.14889   
## Ammonium:StirRate:Temperature  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 ' ' 1

model.final <- aov(Density ~ Ammonium * StirRate, data = powder)
summary(model.final)

##                   Df Sum Sq Mean Sq F value   Pr(>F)    
## Ammonium           1  44.39   44.39   12.17 0.004472 ** 
## StirRate           1  70.69   70.69   19.38 0.000861 ***
## Ammonium:StirRate  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 ' ' 1

The significant effects were Ammonium ((p=0.00447)), Stir Rate ((p=0.000861)), and the Ammonium × Stir Rate interaction ((p=0.01675)).

interaction.plot(powder$StirRate, powder$Ammonium, powder$Density)

2.a.

library(GAD)

density <- c(570,565,583, 1063,1080,1043, 565,510,590,
             528,547,521, 988,1026,1004, 526,538,532)

Position <- rep(c(1,2), each=9)
Temperature <- rep(rep(c(800,825,850), each=3), 2)

Position <- as.fixed(Position)
Temperature <- as.fixed(Temperature)

model.a <- lm(density ~ Position * Temperature)
gad(model.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

The p-values are:

Position: (p=0.00176)

Temperature: (p<0.0001)

Position × Temperature: (p=0.4271)

2.b.

Position <- as.random(Position)
Temperature <- as.random(Temperature)

model.b <- lm(density ~ Position * Temperature)
gad(model.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

The p-values are:

Position: (p=0.05266)

Temperature: (p=0.000865)

Position × Temperature: (p=0.4271)

2.c.

Position <- as.fixed(Position)
Temperature <- as.random(Temperature)

model.c <- lm(density ~ Position * Temperature)
gad(model.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

The p-values are:

Position: (p=0.05266)

Temperature: (p=0.000865)

Position × Temperature: (p=0.4271)

2.d. The p-values differ between the fixed-effects model and the random/mixed-effects models because the appropriate error terms used for the F-tests change. Parts (b) and (c) produced the same p-values. Position was significant only when both factors were treated as fixed, while Temperature was significant in all three models.