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