Problem 1

(a) The model equation is

\[ y_{ijk} = \mu + \alpha _{i} + \beta _{j} +\gamma_{k}+ (\alpha \beta )_{ij}+ (\alpha \gamma )_{ik}+ (\beta \gamma )_{jk}+\alpha\beta\gamma_{ijk} +\epsilon _{ijkl} \]

(b)

library(GAD)
df <- read.csv("https://raw.githubusercontent.com/tmatis12/datafiles/main/PowderProduction.csv")
df$Ammonium <- as.fixed(df$Ammonium)
df$StirRate <- as.fixed(df$StirRate)
df$Temperature <- as.fixed(df$Temperature)

model <- aov(df$Density~df$Ammonium + df$StirRate + df$Temperature + df$Ammonium * df$StirRate + df$Ammonium * df$Temperature + df$StirRate * df$Temperature + df$Ammonium * df$StirRate * df$Temperature)

gad(model)
## $anova
## Analysis of Variance Table
## 
## Response: df$Density
##                                        Df Sum Sq Mean Sq F value   Pr(>F)   
## df$Ammonium                             1 44.389  44.389 11.1803 0.010175 * 
## df$StirRate                             1 70.686  70.686 17.8037 0.002918 **
## df$Temperature                          1  0.328   0.328  0.0826 0.781170   
## df$Ammonium:df$StirRate                 1 28.117  28.117  7.0817 0.028754 * 
## df$Ammonium:df$Temperature              1  0.022   0.022  0.0055 0.942808   
## df$StirRate:df$Temperature              1 10.128  10.128  2.5510 0.148890   
## df$Ammonium:df$StirRate:df$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

Comment: P_value of the three factor interaction is greater than 0.05, which is not significant, and reject the null hypothesis. We dropped the three factor interaction.

model_2<- aov(df$Density~df$Ammonium + df$StirRate + df$Temperature + df$Ammonium * df$StirRate + df$Ammonium * df$Temperature + df$StirRate * df$Temperature)

gad(model_2)
## $anova
## Analysis of Variance Table
## 
## Response: df$Density
##                            Df Sum Sq Mean Sq F value   Pr(>F)   
## df$Ammonium                 1 44.389  44.389 12.0037 0.007109 **
## df$StirRate                 1 70.686  70.686 19.1150 0.001792 **
## df$Temperature              1  0.328   0.328  0.0886 0.772681   
## df$Ammonium:df$StirRate     1 28.117  28.117  7.6033 0.022206 * 
## df$Ammonium:df$Temperature  1  0.022   0.022  0.0059 0.940538   
## df$StirRate:df$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 ' ' 1

Comment: P_value of Ammonium and Temperature is 0.940538, larger than 0.05. we dropped this.

model_3<- aov(df$Density~df$Ammonium + df$StirRate + df$Temperature + df$Ammonium * df$StirRate +  df$StirRate*df$Temperature)

gad(model_3)
## $anova
## Analysis of Variance Table
## 
## Response: df$Density
##                            Df Sum Sq Mean Sq F value    Pr(>F)    
## df$Ammonium                 1 44.389  44.389 13.3287 0.0044560 ** 
## df$StirRate                 1 70.686  70.686 21.2250 0.0009696 ***
## df$Temperature              1  0.328   0.328  0.0984 0.7601850    
## df$Ammonium:df$StirRate     1 28.117  28.117  8.4426 0.0156821 *  
## df$StirRate:df$Temperature  1 10.128  10.128  3.0412 0.1117751    
## Residuals                  10 33.303   3.330                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Comment: P_value of StirRate and Temperature is 0.117, larger than 0.05. we dropped this.

model_4<- aov(df$Density~df$Ammonium + df$StirRate + df$Temperature + df$Ammonium * df$StirRate)

gad(model_4)
## $anova
## Analysis of Variance Table
## 
## Response: df$Density
##                         Df Sum Sq Mean Sq F value   Pr(>F)   
## df$Ammonium              1 44.389  44.389 11.2425 0.006443 **
## df$StirRate              1 70.686  70.686 17.9028 0.001410 **
## df$Temperature           1  0.328   0.328  0.0830 0.778613   
## df$Ammonium:df$StirRate  1 28.117  28.117  7.1211 0.021851 * 
## Residuals               11 43.431   3.948                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Comment: P_value of Temperature is 0.778613, larger than 0.05. we dropped this.

model_5<- aov(df$Density~df$Ammonium + df$StirRate + df$Ammonium * df$StirRate)

gad(model_5)
## $anova
## Analysis of Variance Table
## 
## Response: df$Density
##                         Df Sum Sq Mean Sq F value    Pr(>F)    
## df$Ammonium              1 44.389  44.389 12.1727 0.0044721 ** 
## df$StirRate              1 70.686  70.686 19.3841 0.0008612 ***
## df$Ammonium:df$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

Comment: P_value of every main factors and two way interaction is less than 0.05, all factors now is significant.

The interaction plot

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

Problem 2

Position=c(1,2)
Temperature=c(800,825,850)
df<-expand.grid(Temperature,Position)
df1<-rbind(df,df,df)
colnames(df1)<-c("Temperature", "Position")
df1$Density <- c(
  570, 1063, 565,
  528, 988, 526,
  565, 1080, 510,
  547, 1026, 538, 
  583, 1043, 590,
  521, 1004, 532 )

a)Assume that both Temperature and Position are fixed effects. Report p-values

df1$Temperature <- as.fixed(df1$Temperature)
df1$Position <- as.fixed(df1$Position)
model_6 <- aov(Density~Temperature + Position + Temperature * Position, data = df1)
gad(model_6)
## $anova
## Analysis of Variance Table
## 
## Response: Density
##                      Df Sum Sq Mean Sq  F value   Pr(>F)    
## Temperature           2 945342  472671 1056.117 3.25e-14 ***
## Position              1   7160    7160   15.998 0.001762 ** 
## Temperature:Position  2    818     409    0.914 0.427110    
## Residuals            12   5371     448                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Comment: When both Temperature and Position are fixed effects, the p-value is 0.42711 for two factor interaction, 0.001762 for Position, 3.25e-14 for Temperature.

b) Assume that both Temperature and Position are random effects. Report p-values 

df1$Temperature <- as.random(df1$Temperature)
df1$Position <- as.random(df1$Position)
model_7 <- aov(Density~Temperature + Position + Temperature * Position, data = df1)
gad(model_7)
## $anova
## Analysis of Variance Table
## 
## Response: Density
##                      Df Sum Sq Mean Sq  F value    Pr(>F)    
## Temperature           2 945342  472671 1155.518 0.0008647 ***
## Position              1   7160    7160   17.504 0.0526583 .  
## Temperature:Position  2    818     409    0.914 0.4271101    
## Residuals            12   5371     448                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Comment: When both Temperature and Position are random effects, the p-value is 0.42711 for two factor interaction, 0.0526583 for Position, 0.0008647 for Temperature.

c) Assume the Position effect is fixed and the Temperature effect is random. Report p-values

df1$Temperature <- as.random(df1$Temperature)
df1$Position <- as.fixed(df1$Position)
model_8 <- aov(Density~Temperature + Position + Temperature * Position, data = df1)
gad(model_8)
## $anova
## Analysis of Variance Table
## 
## Response: Density
##                      Df Sum Sq Mean Sq  F value   Pr(>F)    
## Temperature           2 945342  472671 1056.117 3.25e-14 ***
## Position              1   7160    7160   17.504  0.05266 .  
## Temperature:Position  2    818     409    0.914  0.42711    
## Residuals            12   5371     448                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Comment: When Position is fixed effects and Temperature is random effects, the p-value is 0.42711 for two factor interatcion, 0.05266 for Position, 3.25e-14 for Temperature.

d) Comment on similarities and/or differences between the p-values in parts a,b,c.

The P-value for the two factor interaction keeps the same. The p_value for Temperature won’t change no matter is fixed effects or random effects as long as the Position is fixed effects. When they are both random effects, the p value is changed, the temperature is not consistent. The p_value for Position won’t change no matter is fixed effects or random effects as long as the Temperature is random effects.

Complete Code:

library(GAD)
df <- read.csv("https://raw.githubusercontent.com/tmatis12/datafiles/main/PowderProduction.csv")
df$Ammonium <- as.fixed(df$Ammonium)
df$StirRate <- as.fixed(df$StirRate)
df$Temperature <- as.fixed(df$Temperature)

model <- aov(df$Density~df$Ammonium + df$StirRate + df$Temperature + df$Ammonium * df$StirRate + df$Ammonium * df$Temperature + df$StirRate * df$Temperature + df$Ammonium * df$StirRate * df$Temperature)

gad(model)

model_2<- aov(df$Density~df$Ammonium + df$StirRate + df$Temperature + df$Ammonium * df$StirRate + df$Ammonium * df$Temperature + df$StirRate * df$Temperature)

gad(model_2)

model_3<- aov(df$Density~df$Ammonium + df$StirRate + df$Temperature + df$Ammonium * df$StirRate +  df$StirRate*df$Temperature)

gad(model_3)

model_4<- aov(df$Density~df$Ammonium + df$StirRate + df$Temperature + df$Ammonium * df$StirRate)

gad(model_4)

model_5<- aov(df$Density~df$Ammonium + df$StirRate + df$Ammonium * df$StirRate)

gad(model_5)

Position=c(1,2)
Temperature=c(800,825,850)
df<-expand.grid(Temperature,Position)
df1<-rbind(df,df,df)
colnames(df1)<-c("Temperature", "Position")
df1$Density <- c(
  570, 1063, 565,
  528, 988, 526,
  565, 1080, 510,
  547, 1026, 538, 
  583, 1043, 590,
  521, 1004, 532 )

df1$Temperature <- as.fixed(df1$Temperature)
df1$Position <- as.fixed(df1$Position)
model_6 <- aov(Density~Temperature + Position + Temperature * Position, data = df1)
gad(model_6)

df1$Temperature <- as.random(df1$Temperature)
df1$Position <- as.random(df1$Position)
model_7 <- aov(Density~Temperature + Position + Temperature * Position, data = df1)
gad(model_7)

df1$Temperature <- as.random(df1$Temperature)
df1$Position <- as.fixed(df1$Position)
model_8 <- aov(Density~Temperature + Position + Temperature * Position, data = df1)
gad(model_8)