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
\[ 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.
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.
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)