Question 1

library(GAD)
library(DoE.base)
dat <- read.csv("https://raw.githubusercontent.com/tmatis12/datafiles/main/PowderProduction.csv")
dat <- data.frame(dat)
dat
##    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
## 7        30      150           8   12.46
## 8        30      150           8   12.62
## 9         2      100          40   10.95
## 10        2      100          40   17.68
## 11       30      100          40   12.65
## 12       30      100          40   15.96
## 13        2      150          40    8.03
## 14        2      150          40    8.84
## 15       30      150          40   14.96
## 16       30      150          40   14.96

A.The model equation 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} \]

b.

#factors as fixed
dat$Ammonium <- as.fixed(dat$Ammonium)
dat$StirRate <- as.fixed(dat$StirRate)
dat$Temperature <- as.fixed(dat$Temperature)
str(dat)
## 'data.frame':    16 obs. of  4 variables:
##  $ Ammonium   : Factor w/ 2 levels "2","30": 1 1 2 2 1 1 2 2 1 1 ...
##  $ StirRate   : Factor w/ 2 levels "100","150": 1 1 1 1 2 2 2 2 1 1 ...
##  $ Temperature: Factor w/ 2 levels "8","40": 1 1 1 1 1 1 1 1 2 2 ...
##  $ Density    : num  14.68 15.18 15.12 17.48 7.54 ...
model1 <- lm(Density ~ Ammonium*StirRate*Temperature, data=dat)
#coef(model1)
summary(model1)
## 
## Call:
## lm.default(formula = Density ~ Ammonium * StirRate * Temperature, 
##     data = dat)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.3650 -0.4138  0.0000  0.4137  3.3650 
## 
## Coefficients:
##                                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                            14.930      1.409  10.597  5.5e-06 ***
## Ammonium30                              1.370      1.993   0.688  0.51117    
## StirRate150                            -7.830      1.993  -3.930  0.00436 ** 
## Temperature40                          -0.615      1.993  -0.309  0.76547    
## Ammonium30:StirRate150                  4.070      2.818   1.444  0.18664    
## Ammonium30:Temperature40               -1.380      2.818  -0.490  0.63747    
## StirRate150:Temperature40               1.950      2.818   0.692  0.50852    
## Ammonium30:StirRate150:Temperature40    2.465      3.985   0.619  0.55341    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.993 on 8 degrees of freedom
## Multiple R-squared:  0.8301, Adjusted R-squared:  0.6814 
## F-statistic: 5.584 on 7 and 8 DF,  p-value: 0.01359

We eliminate the three order interaction term of Ammonium, Stirrate and Temperature.

model2 <- lm(Density ~ Ammonium + StirRate + Temperature+ Ammonium*StirRate + Ammonium*Temperature +             StirRate*Temperature, data=dat)

summary(model2)
## 
## Call:
## lm.default(formula = Density ~ Ammonium + StirRate + Temperature + 
##     Ammonium * StirRate + Ammonium * Temperature + StirRate * 
##     Temperature, data = dat)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.0569 -0.5969 -0.0950  0.4181  3.6731 
## 
## Coefficients:
##                           Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                15.2381     1.2719  11.980 7.81e-07 ***
## Ammonium30                  0.7537     1.6654   0.453  0.66155    
## StirRate150                -8.4463     1.6654  -5.072  0.00067 ***
## Temperature40              -1.2313     1.6654  -0.739  0.47855    
## Ammonium30:StirRate150      5.3025     1.9230   2.757  0.02221 *  
## Ammonium30:Temperature40   -0.1475     1.9230  -0.077  0.94054    
## StirRate150:Temperature40   3.1825     1.9230   1.655  0.13232    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.923 on 9 degrees of freedom
## Multiple R-squared:  0.822,  Adjusted R-squared:  0.7033 
## F-statistic: 6.926 on 6 and 9 DF,  p-value: 0.005535

We can drop Ammonium and Temperature interaction term.

model3 <- lm(Density ~ Ammonium + StirRate + Temperature+ Ammonium*StirRate + StirRate*Temperature, data=dat)

summary(model3)
## 
## Call:
## lm.default(formula = Density ~ Ammonium + StirRate + Temperature + 
##     Ammonium * StirRate + StirRate * Temperature, data = dat)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.0200 -0.6153 -0.1319  0.3813  3.7100 
## 
## Coefficients:
##                           Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                 15.275      1.117  13.669 8.51e-08 ***
## Ammonium30                   0.680      1.290   0.527 0.609708    
## StirRate150                 -8.446      1.580  -5.344 0.000326 ***
## Temperature40               -1.305      1.290  -1.011 0.335713    
## Ammonium30:StirRate150       5.303      1.825   2.906 0.015682 *  
## StirRate150:Temperature40    3.183      1.825   1.744 0.111775    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.825 on 10 degrees of freedom
## Multiple R-squared:  0.8219, Adjusted R-squared:  0.7328 
## F-statistic: 9.227 on 5 and 10 DF,  p-value: 0.001654

We can drop StirRate and Temperture interaction term looking at the pvalues.

model4 <- lm(Density ~ Ammonium + StirRate + Temperature+ Ammonium*StirRate , data=dat)

summary(model4)
## 
## Call:
## lm.default(formula = Density ~ Ammonium + StirRate + Temperature + 
##     Ammonium * StirRate, data = dat)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.8156 -0.9700  0.1600  0.9638  2.9144 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)             14.4794     1.1108  13.035 4.95e-08 ***
## Ammonium30               0.6800     1.4050   0.484 0.637899    
## StirRate150             -6.8550     1.4050  -4.879 0.000488 ***
## Temperature40            0.2862     0.9935   0.288 0.778613    
## Ammonium30:StirRate150   5.3025     1.9870   2.669 0.021851 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.987 on 11 degrees of freedom
## Multiple R-squared:  0.7677, Adjusted R-squared:  0.6832 
## F-statistic: 9.087 on 4 and 11 DF,  p-value: 0.001703

We drop temperature main effect term.

interaction.plot(dat$StirRate,dat$Ammonium,dat$Density,col=c("BLUE","RED"))

There is interaction between StirRate and Ammonium because the lines are not levelled and they seem to converge.

StirRate and Temperature has significant interaction.

Question 2:

library(agricolae)


A <- c("A1_Organic", "A2_Chemical")     
B <- c("B1_Once",    "B2_Twice")        
C <- c("C1_22C",     "C2_28C")      
library(agricolae)
set.seed(20251106)   
treatments <- expand.grid(Fertilizer = A,
                          Irrigation = B,
                          Temperature = C)


treatments$Label <- with(treatments,
                         paste(Fertilizer, Irrigation, Temperature, sep = "_"))
treatments
##    Fertilizer Irrigation Temperature                       Label
## 1  A1_Organic    B1_Once      C1_22C   A1_Organic_B1_Once_C1_22C
## 2 A2_Chemical    B1_Once      C1_22C  A2_Chemical_B1_Once_C1_22C
## 3  A1_Organic   B2_Twice      C1_22C  A1_Organic_B2_Twice_C1_22C
## 4 A2_Chemical   B2_Twice      C1_22C A2_Chemical_B2_Twice_C1_22C
## 5  A1_Organic    B1_Once      C2_28C   A1_Organic_B1_Once_C2_28C
## 6 A2_Chemical    B1_Once      C2_28C  A2_Chemical_B1_Once_C2_28C
## 7  A1_Organic   B2_Twice      C2_28C  A1_Organic_B2_Twice_C2_28C
## 8 A2_Chemical   B2_Twice      C2_28C A2_Chemical_B2_Twice_C2_28C
design <- design.rcbd(trt   = treatments$Label,
                      r     = 3,        
                      seed  = 20251106,
                      serie = 2)

design$book
##    plots block            treatments$Label
## 1    101     1  A1_Organic_B2_Twice_C2_28C
## 2    102     1 A2_Chemical_B2_Twice_C2_28C
## 3    103     1  A2_Chemical_B1_Once_C2_28C
## 4    104     1   A1_Organic_B1_Once_C1_22C
## 5    105     1  A1_Organic_B2_Twice_C1_22C
## 6    106     1   A1_Organic_B1_Once_C2_28C
## 7    107     1  A2_Chemical_B1_Once_C1_22C
## 8    108     1 A2_Chemical_B2_Twice_C1_22C
## 9    201     2 A2_Chemical_B2_Twice_C1_22C
## 10   202     2  A2_Chemical_B1_Once_C1_22C
## 11   203     2   A1_Organic_B1_Once_C2_28C
## 12   204     2   A1_Organic_B1_Once_C1_22C
## 13   205     2  A1_Organic_B2_Twice_C2_28C
## 14   206     2  A1_Organic_B2_Twice_C1_22C
## 15   207     2  A2_Chemical_B1_Once_C2_28C
## 16   208     2 A2_Chemical_B2_Twice_C2_28C
## 17   301     3  A2_Chemical_B1_Once_C2_28C
## 18   302     3 A2_Chemical_B2_Twice_C2_28C
## 19   303     3  A1_Organic_B2_Twice_C1_22C
## 20   304     3 A2_Chemical_B2_Twice_C1_22C
## 21   305     3  A2_Chemical_B1_Once_C1_22C
## 22   306     3  A1_Organic_B2_Twice_C2_28C
## 23   307     3   A1_Organic_B1_Once_C1_22C
## 24   308     3   A1_Organic_B1_Once_C2_28C
library(dplyr)
design$book %>%
  arrange(block, plots)
##    plots block            treatments$Label
## 1    101     1  A1_Organic_B2_Twice_C2_28C
## 2    102     1 A2_Chemical_B2_Twice_C2_28C
## 3    103     1  A2_Chemical_B1_Once_C2_28C
## 4    104     1   A1_Organic_B1_Once_C1_22C
## 5    105     1  A1_Organic_B2_Twice_C1_22C
## 6    106     1   A1_Organic_B1_Once_C2_28C
## 7    107     1  A2_Chemical_B1_Once_C1_22C
## 8    108     1 A2_Chemical_B2_Twice_C1_22C
## 9    201     2 A2_Chemical_B2_Twice_C1_22C
## 10   202     2  A2_Chemical_B1_Once_C1_22C
## 11   203     2   A1_Organic_B1_Once_C2_28C
## 12   204     2   A1_Organic_B1_Once_C1_22C
## 13   205     2  A1_Organic_B2_Twice_C2_28C
## 14   206     2  A1_Organic_B2_Twice_C1_22C
## 15   207     2  A2_Chemical_B1_Once_C2_28C
## 16   208     2 A2_Chemical_B2_Twice_C2_28C
## 17   301     3  A2_Chemical_B1_Once_C2_28C
## 18   302     3 A2_Chemical_B2_Twice_C2_28C
## 19   303     3  A1_Organic_B2_Twice_C1_22C
## 20   304     3 A2_Chemical_B2_Twice_C1_22C
## 21   305     3  A2_Chemical_B1_Once_C1_22C
## 22   306     3  A1_Organic_B2_Twice_C2_28C
## 23   307     3   A1_Organic_B1_Once_C1_22C
## 24   308     3   A1_Organic_B1_Once_C2_28C
write.csv(design$book, "TomatoYield_DesignBook.csv", row.names = FALSE)