knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning = FALSE)

Problem 1

1a) Model

model: \[ y = \mu + A + B + C + AB + AC + BC + ABC + \varepsilon \]

1b) Analysis in R

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

# Treat as factors 
powder$Ammonium    <- as.factor(powder$Ammonium)
powder$StirRate    <- as.factor(powder$StirRate)
powder$Temperature <- as.factor(powder$Temperature)

# Full model
fit_full <- aov(Density ~ Ammonium*StirRate*Temperature, data = powder)
anova(fit_full)
## 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
# Step down reduction 
fit1 <- update(fit_full, . ~ . - Ammonium:StirRate:Temperature)
drop1(fit1, test="F")
## Single term deletions
## 
## Model:
## Density ~ Ammonium + StirRate + Temperature + Ammonium:StirRate + 
##     Ammonium:Temperature + StirRate:Temperature
##                      Df Sum of Sq    RSS    AIC F value  Pr(>F)  
## <none>                            33.281 25.719                  
## Ammonium:StirRate     1   28.1165 61.398 33.517  7.6033 0.02221 *
## Ammonium:Temperature  1    0.0218 33.303 23.729  0.0059 0.94054  
## StirRate:Temperature  1   10.1283 43.410 27.970  2.7389 0.13232  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
fit2 <- update(fit1, . ~ . - Ammonium:Temperature)
drop1(fit2, test="F")
## Single term deletions
## 
## Model:
## Density ~ Ammonium + StirRate + Temperature + Ammonium:StirRate + 
##     StirRate:Temperature
##                      Df Sum of Sq    RSS    AIC F value  Pr(>F)  
## <none>                            33.303 23.729                  
## Ammonium:StirRate     1    28.116 61.420 31.522  8.4426 0.01568 *
## StirRate:Temperature  1    10.128 43.431 25.977  3.0412 0.11178  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
fit3 <- update(fit2, . ~ . - StirRate:Temperature)
drop1(fit3, test="F")
## Single term deletions
## 
## Model:
## Density ~ Ammonium + StirRate + Temperature + Ammonium:StirRate
##                   Df Sum of Sq    RSS    AIC F value  Pr(>F)  
## <none>                         43.431 25.978                  
## Temperature        1    0.3278 43.759 24.098  0.0830 0.77861  
## Ammonium:StirRate  1   28.1165 71.548 31.964  7.1211 0.02185 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# what stays significant when run: A, B, and A:B
fit_final <- aov(Density ~ Ammonium*StirRate, data = powder)
anova(fit_final)
## 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
summary(fit_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
# Quick interaction plots 
par(mfrow=c(1,3))
with(powder, interaction.plot(Ammonium, StirRate, Density, main="A:B"))
with(powder, interaction.plot(Ammonium, Temperature, Density, main="A:C"))
with(powder, interaction.plot(StirRate, Temperature, Density, main="B:C"))

par(mfrow=c(1,1))

# Quick residual checks
par(mfrow=c(2,2)); plot(fit_final); par(mfrow=c(1,1))

Problem 2 — 2^3 RCBD

Randomize 8 treatments within each of 3 blocks.

set.seed(271828)

A <- c("A1","A2")
B <- c("B1","B2")
C <- c("C1","C2")

trts <- with(expand.grid(A=A,B=B,C=C), paste(A,B,C,sep=":"))

make_block <- function(bname, trts) {
  data.frame(block=bname,
             plot=1:length(trts),
             trt=sample(trts))
}

book <- rbind(
  make_block("Block1", trts),
  make_block("Block2", trts),
  make_block("Block3", trts)
)

book
##     block plot      trt
## 1  Block1    1 A2:B1:C2
## 2  Block1    2 A1:B2:C2
## 3  Block1    3 A2:B2:C1
## 4  Block1    4 A2:B2:C2
## 5  Block1    5 A1:B1:C1
## 6  Block1    6 A2:B1:C1
## 7  Block1    7 A1:B2:C1
## 8  Block1    8 A1:B1:C2
## 9  Block2    1 A1:B1:C2
## 10 Block2    2 A2:B2:C1
## 11 Block2    3 A1:B2:C1
## 12 Block2    4 A1:B1:C1
## 13 Block2    5 A2:B1:C1
## 14 Block2    6 A1:B2:C2
## 15 Block2    7 A2:B1:C2
## 16 Block2    8 A2:B2:C2
## 17 Block3    1 A1:B2:C2
## 18 Block3    2 A2:B1:C2
## 19 Block3    3 A1:B1:C2
## 20 Block3    4 A1:B2:C1
## 21 Block3    5 A2:B2:C1
## 22 Block3    6 A2:B1:C1
## 23 Block3    7 A2:B2:C2
## 24 Block3    8 A1:B1:C1