library(GAD)
## Loading required package: matrixStats
## Loading required package: R.methodsS3
## R.methodsS3 v1.8.2 (2022-06-13 22:00:14 UTC) successfully loaded. See ?R.methodsS3 for help.
dat <- read.csv("https://raw.githubusercontent.com/tmatis12/datafiles/main/PowderProduction.csv")
The equation for a three factor analysis is the following:
\[ y_{ijlk}=\mu+\alpha_i+\beta_j+\gamma_l+\alpha\beta_{ij}+\alpha\gamma_{il}+\beta\gamma_{jl}+\alpha\beta\gamma_{ijl}+\epsilon_{ijlk} \]
The hypothesis that we are trying to reject are the following (not in order of relevance):
\[ Ho: \alpha\beta\gamma_{ijl}= 0\\ Ha: \alpha\beta\gamma_{ijl}\neq 0 \]
\[ Ho: \alpha\beta_{ij}= 0\\ Ha: \alpha\beta_{ij}\neq 0 \]
\[ Ho: \beta\gamma_{jl}= 0\\ Ha: \beta\gamma_{jl}\neq 0 \]
\[ Ho: \alpha\gamma_{il}= 0\\ Ha: \alpha\gamma_{il}\neq 0 \]
\[ Ho: \alpha_{i}= 0\\ Ha: \alpha_{i}\neq 0 \]
\[ Ho: \beta_{j}= 0\\ Ha: \beta_{j}\neq 0 \]
\[ Ho: \gamma_{l}= 0 \\ Ha: \gamma_{l}\neq 0 \]
dat$Ammonium <- as.fixed(dat$Ammonium)
dat$StirRate <- as.fixed(dat$StirRate)
dat$Temperature <- as.fixed(dat$Temperature)
To test our hypothesis, let’s perform the following:
\[ Ho: \alpha\beta\gamma_{ijl}= 0\\ Ha: \alpha\beta\gamma_{ijl}\neq 0 \]
model.aov <- aov(Density~Ammonium+StirRate+Temperature+Ammonium*StirRate+Ammonium*Temperature+StirRate*Temperature+Ammonium*StirRate*Temperature, data=dat)
summary(model.aov)
## Df Sum Sq Mean Sq F value Pr(>F)
## Ammonium 1 44.39 44.39 11.180 0.01018 *
## StirRate 1 70.69 70.69 17.804 0.00292 **
## Temperature 1 0.33 0.33 0.083 0.78117
## Ammonium:StirRate 1 28.12 28.12 7.082 0.02875 *
## Ammonium:Temperature 1 0.02 0.02 0.005 0.94281
## StirRate:Temperature 1 10.13 10.13 2.551 0.14889
## Ammonium:StirRate:Temperature 1 1.52 1.52 0.383 0.55341
## Residuals 8 31.76 3.97
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
From the results above, it we cannot reject Ho, since the p-value is greater than \(\alpha\). Therefore, the interaction between the three factors is not significant.
Now, let’s test the least significant effect in order to validate the two factors interaction by removing the three factor interaction.
\[ Ho: \alpha\gamma_{il}= 0\\ Ha: \alpha\gamma_{il}\neq 0 \]
model.aov <- aov(Density~Ammonium+StirRate+Temperature+Ammonium*StirRate+Ammonium*Temperature+StirRate*Temperature, data=dat)
summary(model.aov)
## Df Sum Sq Mean Sq F value Pr(>F)
## Ammonium 1 44.39 44.39 12.004 0.00711 **
## StirRate 1 70.69 70.69 19.115 0.00179 **
## Temperature 1 0.33 0.33 0.089 0.77268
## Ammonium:StirRate 1 28.12 28.12 7.603 0.02221 *
## Ammonium:Temperature 1 0.02 0.02 0.006 0.94054
## StirRate:Temperature 1 10.13 10.13 2.739 0.13232
## Residuals 9 33.28 3.70
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Since the p-value is greater than \(\alpha\), we cannot reject Ho. Therefore, we need to check the next least significant interaction by removing alpha gamma interation.
\[ Ho: \beta\gamma_{jl}= 0\\ Ha: \beta\gamma_{jl}\neq 0 \]
model.aov <- aov(Density~Ammonium+StirRate+Temperature+Ammonium*StirRate+StirRate*Temperature, data=dat)
summary(model.aov)
## Df Sum Sq Mean Sq F value Pr(>F)
## Ammonium 1 44.39 44.39 13.329 0.00446 **
## StirRate 1 70.69 70.69 21.225 0.00097 ***
## Temperature 1 0.33 0.33 0.098 0.76019
## Ammonium:StirRate 1 28.12 28.12 8.443 0.01568 *
## StirRate:Temperature 1 10.13 10.13 3.041 0.11178
## Residuals 10 33.30 3.33
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Since the p-value is greater than \(\alpha\), we cannot reject Ho. Therefore, we need to test the last two factor interaction to see its relevance by removing beta gamma interaction.
\[ Ho: \alpha\beta_{ij}= 0\\ Ha: \alpha\beta_{ij}\neq 0 \]
model.aov <- aov(Density~Ammonium+StirRate+Temperature+Ammonium*StirRate, data=dat)
summary(model.aov)
## Df Sum Sq Mean Sq F value Pr(>F)
## Ammonium 1 44.39 44.39 11.242 0.00644 **
## StirRate 1 70.69 70.69 17.903 0.00141 **
## Temperature 1 0.33 0.33 0.083 0.77861
## Ammonium:StirRate 1 28.12 28.12 7.121 0.02185 *
## Residuals 11 43.43 3.95
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Since the p-value of the last two factor interaction is lesser than \(\alpha\) (\(p-value=0.022<\alpha=0.05\)), we can reject Ho and consider that there is a significant interaction between the Ammonium percentage and the Stir Rate.
Since there is an interaction going on, we can check the plot by
interaction.plot(dat$Ammonium,dat$StirRate,dat$Density)
From the plot, since the curves are not parallel, we can confirm the interaction between the two factors.
\[ Ho: \alpha\beta_{ij}=0 \\ Ha: \alpha\beta_{ij}\neq0 \]
library(GAD)
position <- c(rep(1,9),rep(2,9))
temperature <- c(rep(seq(800,850,25),6))
response <- c(570, 1063, 565,
565, 1080, 510,
583, 1043, 590,
528, 988, 526,
547, 1026, 538,
521, 1004, 532)
position <- as.fixed(position)
temperature <- as.fixed(temperature)
dat <- data.frame(position,temperature,response)
dat
## position temperature response
## 1 1 800 570
## 2 1 825 1063
## 3 1 850 565
## 4 1 800 565
## 5 1 825 1080
## 6 1 850 510
## 7 1 800 583
## 8 1 825 1043
## 9 1 850 590
## 10 2 800 528
## 11 2 825 988
## 12 2 850 526
## 13 2 800 547
## 14 2 825 1026
## 15 2 850 538
## 16 2 800 521
## 17 2 825 1004
## 18 2 850 532
model1 <- aov(response~position+temperature+position*temperature,data=dat)
gad(model1)
## Analysis of Variance Table
##
## Response: response
## Df Sum Sq Mean Sq F value Pr(>F)
## position 1 7160 7160 15.998 0.001762 **
## temperature 2 945342 472671 1056.117 3.25e-14 ***
## position:temperature 2 818 409 0.914 0.427110
## Residual 12 5371 448
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
P value of the interaction is 0.427110 which is larger than \(\alpha = 0.05\) level of significance, so we don’t reject H0 hence there is no two factors interaction present. Also p value for the position and temperature is 0.001762 and 3.25e-14 respectively.
\[ Ho: \alpha\beta_{ij}=0 \\ Ha: \alpha\beta_{ij}\neq0 \]
library(GAD)
position <- c(rep(1,9),rep(2,9))
temperature <- c(rep(seq(800,850,25),6))
response <- c(570, 1063, 565,
565, 1080, 510,
583, 1043, 590,
528, 988, 526,
547, 1026, 538,
521, 1004, 532)
position <- as.random(position)
temperature <- as.random(temperature)
dat2 <- data.frame(position,temperature,response)
dat2
## position temperature response
## 1 1 800 570
## 2 1 825 1063
## 3 1 850 565
## 4 1 800 565
## 5 1 825 1080
## 6 1 850 510
## 7 1 800 583
## 8 1 825 1043
## 9 1 850 590
## 10 2 800 528
## 11 2 825 988
## 12 2 850 526
## 13 2 800 547
## 14 2 825 1026
## 15 2 850 538
## 16 2 800 521
## 17 2 825 1004
## 18 2 850 532
model2 <- aov(response~position+temperature+position*temperature,data=dat2)
gad(model2)
## Analysis of Variance Table
##
## Response: response
## Df Sum Sq Mean Sq F value Pr(>F)
## position 1 7160 7160 17.504 0.0526583 .
## temperature 2 945342 472671 1155.518 0.0008647 ***
## position:temperature 2 818 409 0.914 0.4271101
## Residual 12 5371 448
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
P value of the interaction is 0.4271101 which is larger than \(\alpha = 0.05\) level of significance, so we don’t reject H0 hence there is no two factors interaction present. Also p value for the position and temperature is 0.0526583 and 0.0008647 respectively.
\[ Ho: \alpha\beta_{ij}=0 \\ Ha: \alpha\beta_{ij}\neq0 \]
library(GAD)
position <- c(rep(1,9),rep(2,9))
temperature <- c(rep(seq(800,850,25),6))
response <- c(570, 1063, 565,
565, 1080, 510,
583, 1043, 590,
528, 988, 526,
547, 1026, 538,
521, 1004, 532)
position <- as.fixed(position)
temperature <- as.random(temperature)
dat3 <- data.frame(position,temperature,response)
dat3
## position temperature response
## 1 1 800 570
## 2 1 825 1063
## 3 1 850 565
## 4 1 800 565
## 5 1 825 1080
## 6 1 850 510
## 7 1 800 583
## 8 1 825 1043
## 9 1 850 590
## 10 2 800 528
## 11 2 825 988
## 12 2 850 526
## 13 2 800 547
## 14 2 825 1026
## 15 2 850 538
## 16 2 800 521
## 17 2 825 1004
## 18 2 850 532
model3 <- aov(response~position+temperature+position*temperature,data=dat3)
gad(model3)
## Analysis of Variance Table
##
## Response: response
## Df Sum Sq Mean Sq F value Pr(>F)
## position 1 7160 7160 17.504 0.05266 .
## temperature 2 945342 472671 1056.117 3.25e-14 ***
## position:temperature 2 818 409 0.914 0.42711
## Residual 12 5371 448
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
P value of the interaction is 0.427110 which is larger than \(\alpha = 0.05\) level of significance, so we don’t reject H0 hence there is no two factors interaction present. Also p value for the position and temperature is 0.05266 and 3.25e-14 respectively.
Comments on p values:
The interaction term is not significant for all the part A, B and C.
When both factors are fixed, position and temperature both are significant at alpha= 0.01 and 0.001 level of significance.
When both factors are random, position and temperature both are significant at alpha= 0.1 and 0.001 level of significance.
When the factors have mixed effect, position and temperature both are significant at alpha= 0.1 and 0.001 level of significance.
But, it is possible to see a difference in the p-values of the main effects:
For fixed models, we are getting different p-values when compared with the random models. However, for mixed models, we are getting smaller p-values for the main effects.
Temperature p-value is the same for Fixed and Mixed Models, but different in the random
Position p-value is the same for mixed and Random but different for fixed