Question 1

library(GAD)
library(DoE.base)
## Warning: package 'DoE.base' was built under R version 4.5.2
## Loading required package: grid
## Loading required package: conf.design
## Warning: package 'conf.design' was built under R version 4.5.2
## Registered S3 method overwritten by 'DoE.base':
##   method           from       
##   factorize.factor conf.design
## 
## Attaching package: 'DoE.base'
## The following objects are masked from 'package:stats':
## 
##     aov, lm
## The following object is masked from 'package:graphics':
## 
##     plot.design
## The following object is masked from 'package:base':
## 
##     lengths
# loading the data 
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

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} \]

#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 because the p-value is more than the alpha(0.05).

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 because it is not significant(>alpha).

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 because it is not signinficant.

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 because it is not significant.

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 Ammonium has significant interaction.

Question 2:

library(agricolae)

?design.ab
## starting httpd help server ... done
factorA <- c("A_1 Organic compost", "A_2 Chemical fertilizer")
factorB <- c("B_2 Once per day", "B_2 Twice per day")
factorC <- c("C_1 22°C", "C_2 28°C")

# factorial treatment structure
trts <- list(A = factorA, B = factorB, C = factorC)
trts<- c(2,2,2)

# Generate a 2^3 factorial design with 3 blocks (replications)
design <- design.ab(trt=trts, r = 3, design = "rcbd", seed = 20255)

design$book
##    plots block A B C
## 1    101     1 1 1 2
## 2    102     1 1 2 2
## 3    103     1 2 2 1
## 4    104     1 1 2 1
## 5    105     1 2 2 2
## 6    106     1 2 1 2
## 7    107     1 1 1 1
## 8    108     1 2 1 1
## 9    109     2 1 1 2
## 10   110     2 1 2 2
## 11   111     2 2 2 2
## 12   112     2 2 1 2
## 13   113     2 2 2 1
## 14   114     2 2 1 1
## 15   115     2 1 2 1
## 16   116     2 1 1 1
## 17   117     3 1 1 2
## 18   118     3 2 2 1
## 19   119     3 2 2 2
## 20   120     3 1 2 2
## 21   121     3 2 1 2
## 22   122     3 1 1 1
## 23   123     3 1 2 1
## 24   124     3 2 1 1
colnames(design$book)<- c("plot","block","Fertilizer Type","Irrigation Frequency", "Temperature Setting")
design$book
##    plot block Fertilizer Type Irrigation Frequency Temperature Setting
## 1   101     1               1                    1                   2
## 2   102     1               1                    2                   2
## 3   103     1               2                    2                   1
## 4   104     1               1                    2                   1
## 5   105     1               2                    2                   2
## 6   106     1               2                    1                   2
## 7   107     1               1                    1                   1
## 8   108     1               2                    1                   1
## 9   109     2               1                    1                   2
## 10  110     2               1                    2                   2
## 11  111     2               2                    2                   2
## 12  112     2               2                    1                   2
## 13  113     2               2                    2                   1
## 14  114     2               2                    1                   1
## 15  115     2               1                    2                   1
## 16  116     2               1                    1                   1
## 17  117     3               1                    1                   2
## 18  118     3               2                    2                   1
## 19  119     3               2                    2                   2
## 20  120     3               1                    2                   2
## 21  121     3               2                    1                   2
## 22  122     3               1                    1                   1
## 23  123     3               1                    2                   1
## 24  124     3               2                    1                   1