Load packages

library(readxl)
library(ExpDes)
library(agricolae)
library(tidyverse)
library(ggstatsplot)
library(car)
library(patchwork)
library(emmeans)
library(multcomp)

Three-factor factorial experiment

three_factor <- read_excel("three_factorial.xlsx")

three_factor <- three_factor %>% 
  mutate(Block = factor(Block),
         A = factor(A, labels = c("A1", "A2")),
         B = factor(B, labels = c("B1", "B2")),
         C = factor(C, labels=c("C1", "C2", "C3")))
with(three_factor, fat3.rbd(factor1 = A,
                            factor2 = B, 
                            factor3 = C, 
                            block = Block, 
                            resp = Y, 
                            quali = c(TRUE, TRUE, TRUE),
                            mcomp = "tukey", 
                            fac.names = c("A", "B", "C")))
## ------------------------------------------------------------------------
## Legend:
## FACTOR 1:  A 
## FACTOR 2:  B 
## FACTOR 3:  C 
## ------------------------------------------------------------------------
## 
## ------------------------------------------------------------------------
## Analysis of Variance Table
## ------------------------------------------------------------------------
##           DF         SS         MS      Fc  Pr>Fc
## Block      5  131.89347   26.37869  1.4099 0.2351
## A          1  313.35877  313.35877  16.749  1e-04
## B          1 1793.55576 1793.55576 95.8655      0
## C          2  283.16803  141.58402  7.5677 0.0012
## A*B        1    2.73058    2.73058  0.1459 0.7039
## A*C        2  129.18232   64.59116  3.4524 0.0387
## B*C        2   27.59502   13.79751  0.7375  0.483
## A*B*C      2  286.59562  143.29781  7.6593 0.0012
## Residuals 55 1028.99941   18.70908               
## Total     66 3997.07898                          
## ------------------------------------------------------------------------
## CV = 5.57 %
## 
## ------------------------------------------------------------------------
## Shapiro-Wilk normality test
## p-value:  0.1545518 
## According to Shapiro-Wilk normality test at 5% of significance, residuals can be considered normal.
## ------------------------------------------------------------------------
## 
## 
## 
## Significant A*B*C  interaction: analyzing the interaction
## ------------------------------------------------------------------------
## 
## Analyzing  A  inside of each level of  B and C 
## ------------------------------------------------------------------------
## ------------------------------------------------------------------------
## Analysis of Variance Table
## ------------------------------------------------------------------------
##           DF          SS         MS        Fc    Pr>Fc
## A: B1 C1   1   10.779355  10.779355  0.576156 0.451064
## A: B1 C2   1   68.454804  68.454804  3.658908 0.060982
## A: B1 C3   1   65.603549  65.603549  3.506509 0.066448
## A: B2 C1   1  577.446016 577.446016 30.864479    1e-06
## A: B2 C2   1    4.080335   4.080335  0.218094 0.642341
## A: B2 C3   1    5.503246   5.503246  0.294148 0.589765
## Residuals 55 1028.999406  18.709080                   
## ------------------------------------------------------------------------
## 
## 
## 
##  A  inside of the combination of the levels  B1  of  B  and  C1  of  C 
## 
## According to the F test, the means of this factor are statistical equal.
## ------------------------------------------------------------------------
##     Levels     Means
## 1       A1  76.05178
## 2       A2  74.15623
## ------------------------------------------------------------------------
## 
## 
##  A  inside of the combination of the levels  B1  of  B  and  C2  of  C 
## 
## According to the F test, the means of this factor are statistical equal.
## ------------------------------------------------------------------------
##     Levels     Means
## 1       A1  73.13877
## 2       A2  68.36193
## ------------------------------------------------------------------------
## 
## 
##  A  inside of the combination of the levels  B1  of  B  and  C3  of  C 
## 
## According to the F test, the means of this factor are statistical equal.
## ------------------------------------------------------------------------
##     Levels     Means
## 1       A1  74.45586
## 2       A2  69.77956
## ------------------------------------------------------------------------
## 
## 
##  A  inside of the combination of the levels  B2  of  B  and  C1  of  C 
## ------------------------------------------------------------------------
## Tukey's test
## ------------------------------------------------------------------------
## Groups Treatments Means
## a     A1      92.73885 
##  b    A2      78.86506 
## ------------------------------------------------------------------------
## 
## 
##  A  inside of the combination of the levels  B2  of  B  and  C2  of  C 
## 
## According to the F test, the means of this factor are statistical equal.
## ------------------------------------------------------------------------
##     Levels     Means
## 1       A1  82.34154
## 2       A2  81.17530
## ------------------------------------------------------------------------
## 
## 
##  A  inside of the combination of the levels  B2  of  B  and  C3  of  C 
## 
## According to the F test, the means of this factor are statistical equal.
## ------------------------------------------------------------------------
##     Levels     Means
## 1       A1  79.68074
## 2       A2  81.03514
## ------------------------------------------------------------------------
## 
## 
## 
## Analyzing  B  inside of each level of  A and C 
## ------------------------------------------------------------------------
## ------------------------------------------------------------------------
## Analysis of Variance Table
## ------------------------------------------------------------------------
##           DF         SS        MS        Fc    Pr>Fc
## B: A1 C1   1  835.37433 835.37433 44.650743        0
## B: A1 C2   1  254.07290 254.07290 13.580192 0.000524
## B: A1 C3   1   81.89790  81.89790  4.377441 0.041046
## B: A2 C1   1   66.51918  66.51918  3.555449 0.064636
## B: A2 C2   1  492.54799 492.54799 26.326681    4e-06
## B: A2 C3   1  380.06468 380.06468 20.314451  3.5e-05
## Residuals 55 1028.99941  18.70908                   
## ------------------------------------------------------------------------
## 
## 
## 
##  B  inside of the combination of the levels  A1  of  A  and  C1  of  C 
## ------------------------------------------------------------------------
## Tukey's test
## ------------------------------------------------------------------------
## Groups Treatments Means
## a     B2      92.73885 
##  b    B1      76.05178 
## ------------------------------------------------------------------------
## 
## 
##  B  inside of the combination of the levels  A1  of  A  and  C2  of  C 
## ------------------------------------------------------------------------
## Tukey's test
## ------------------------------------------------------------------------
## Groups Treatments Means
## a     B2      82.34154 
##  b    B1      73.13877 
## ------------------------------------------------------------------------
## 
## 
##  B  inside of the combination of the levels  A1  of  A  and  C3  of  C 
## ------------------------------------------------------------------------
## Tukey's test
## ------------------------------------------------------------------------
## Groups Treatments Means
## a     B2      79.68074 
##  b    B1      74.45586 
## ------------------------------------------------------------------------
## 
## 
##  B  inside of the combination of the levels  A2  of  A  and  C1  of  C 
## 
## According to the F test, the means of this factor are statistical equal.
## ------------------------------------------------------------------------
##     Levels     Means
## 1       B1  74.15623
## 2       B2  78.86506
## ------------------------------------------------------------------------
## 
## 
##  B  inside of the combination of the levels  A2  of  A  and  C2  of  C 
## ------------------------------------------------------------------------
## Tukey's test
## ------------------------------------------------------------------------
## Groups Treatments Means
## a     B2      81.1753 
##  b    B1      68.36193 
## ------------------------------------------------------------------------
## 
## 
##  B  inside of the combination of the levels  A2  of  A  and  C3  of  C 
## ------------------------------------------------------------------------
## Tukey's test
## ------------------------------------------------------------------------
## Groups Treatments Means
## a     B2      81.03514 
##  b    B1      69.77956 
## ------------------------------------------------------------------------
## 
## Analyzing  C  inside of each level of  A and B 
## ------------------------------------------------------------------------
## ------------------------------------------------------------------------
## Analysis of Variance Table
## ------------------------------------------------------------------------
##           DF         SS        MS        Fc    Pr>Fc
## C: A1 B1   2   25.53464  12.76732  0.682413 0.509623
## C: A1 B2   2  571.39621 285.69811 15.270559    5e-06
## C: A2 B1   2  109.47787  54.73893  2.925795 0.062015
## C: A2 B2   2   20.13228  10.06614  0.538035 0.586936
## Residuals 55 1028.99941  18.70908                   
## ------------------------------------------------------------------------
## 
## 
## 
##  C  inside of the combination of the levels  A1  of  A  and  B1  of  B 
## 
## According to the F test, the means of this factor are statistical equal.
## ------------------------------------------------------------------------
##     Levels     Means
## 1       C1  76.05178
## 2       C2  73.13877
## 3       C3  74.45586
## ------------------------------------------------------------------------
## 
## 
##  C  inside of the combination of the levels  A1  of  A  and  B2  of  B 
## ------------------------------------------------------------------------
## Tukey's test
## ------------------------------------------------------------------------
## Groups Treatments Means
## a     C1      92.73885 
##  b    C2      82.34154 
##  b    C3      79.68074 
## ------------------------------------------------------------------------
## 
## 
##  C  inside of the combination of the levels  A2  of  A  and  B1  of  B 
## 
## According to the F test, the means of this factor are statistical equal.
## ------------------------------------------------------------------------
##     Levels     Means
## 1       C1  74.15623
## 2       C2  68.36193
## 3       C3  69.77956
## ------------------------------------------------------------------------
## 
## 
##  C  inside of the combination of the levels  A2  of  A  and  B2  of  B 
## 
## According to the F test, the means of this factor are statistical equal.
## ------------------------------------------------------------------------
##     Levels     Means
## 1       C1  78.86506
## 2       C2  81.17530
## 3       C3  81.03514
## ------------------------------------------------------------------------

ANOVA Table using aov()

threefactor.model <- aov(Y ~ Block + A*B*C,
                         data = three_factor)

anova(threefactor.model)
## Analysis of Variance Table
## 
## Response: Y
##           Df  Sum Sq Mean Sq F value    Pr(>F)    
## Block      5  131.89   26.38  1.4099 0.2351306    
## A          1  313.36  313.36 16.7490 0.0001408 ***
## B          1 1793.56 1793.56 95.8655 1.182e-13 ***
## C          2  283.17  141.58  7.5677 0.0012494 ** 
## A:B        1    2.73    2.73  0.1459 0.7039086    
## A:C        2  129.18   64.59  3.4524 0.0386862 *  
## B:C        2   27.60   13.80  0.7375 0.4829890    
## A:B:C      2  286.60  143.30  7.6593 0.0011629 ** 
## Residuals 55 1029.00   18.71                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Compute 3-way treatment combination means and assign letters

abc.means <- three_factor  %>%  
  group_by(A,B,C)  %>%  
  dplyr::summarise(Mean = mean(Y),
            SD = sd(Y),
            n = length(Y))  %>%  
  dplyr::mutate(SE = SD/sqrt(n))

abc.means
## # A tibble: 12 × 7
## # Groups:   A, B [4]
##    A     B     C      Mean    SD     n    SE
##    <fct> <fct> <fct> <dbl> <dbl> <int> <dbl>
##  1 A1    B1    C1     76.1  3.85     6  1.57
##  2 A1    B1    C2     73.1  4.77     6  1.95
##  3 A1    B1    C3     74.5  4.89     6  2.00
##  4 A1    B2    C1     92.7  5.12     6  2.09
##  5 A1    B2    C2     82.3  5.00     6  2.04
##  6 A1    B2    C3     79.7  4.05     6  1.65
##  7 A2    B1    C1     74.2  3.69     6  1.51
##  8 A2    B1    C2     68.4  4.08     6  1.67
##  9 A2    B1    C3     69.8  2.72     6  1.11
## 10 A2    B2    C1     78.9  5.32     6  2.17
## 11 A2    B2    C2     81.2  4.62     6  1.89
## 12 A2    B2    C3     81.0  3.98     6  1.63
ABC.means <- emmeans(threefactor.model, specs=c("A", "B", "C"))
ABC.means
##  A  B  C  emmean   SE df lower.CL upper.CL
##  A1 B1 C1   76.1 1.77 55     72.5     79.6
##  A2 B1 C1   74.2 1.77 55     70.6     77.7
##  A1 B2 C1   92.7 1.77 55     89.2     96.3
##  A2 B2 C1   78.9 1.77 55     75.3     82.4
##  A1 B1 C2   73.1 1.77 55     69.6     76.7
##  A2 B1 C2   68.4 1.77 55     64.8     71.9
##  A1 B2 C2   82.3 1.77 55     78.8     85.9
##  A2 B2 C2   81.2 1.77 55     77.6     84.7
##  A1 B1 C3   74.5 1.77 55     70.9     78.0
##  A2 B1 C3   69.8 1.77 55     66.2     73.3
##  A1 B2 C3   79.7 1.77 55     76.1     83.2
##  A2 B2 C3   81.0 1.77 55     77.5     84.6
## 
## Results are averaged over the levels of: Block 
## Confidence level used: 0.95
#AGenerate the letters
ABC.means.cld <- cld(ABC.means,
                   Letters=letters,
                   decreasing = TRUE)

ABC.means.cld <- ABC.means.cld  %>%  
  arrange(A, B, C)  %>%  
  mutate(SE1 = abc.means$SE)

ABC.means.cld
##     A  B  C   emmean       SE df lower.CL upper.CL .group      SE1
## 1  A1 B1 C1 76.05178 1.765837 55 72.51297 79.59060    bcd 1.573747
## 2  A1 B1 C2 73.13877 1.765837 55 69.59996 76.67759     cd 1.945328
## 3  A1 B1 C3 74.45586 1.765837 55 70.91705 77.99468    bcd 1.995871
## 4  A1 B2 C1 92.73885 1.765837 55 89.20003 96.27766   a    2.088637
## 5  A1 B2 C2 82.34154 1.765837 55 78.80273 85.88036    b   2.041608
## 6  A1 B2 C3 79.68074 1.765837 55 76.14192 83.21955    bc  1.651725
## 7  A2 B1 C1 74.15623 1.765837 55 70.61741 77.69505    bcd 1.506547
## 8  A2 B1 C2 68.36193 1.765837 55 64.82311 71.90074      d 1.666405
## 9  A2 B1 C3 69.77956 1.765837 55 66.24074 73.31837      d 1.110290
## 10 A2 B2 C1 78.86506 1.765837 55 75.32624 82.40387    bc  2.170447
## 11 A2 B2 C2 81.17530 1.765837 55 77.63649 84.71412    bc  1.885857
## 12 A2 B2 C3 81.03514 1.765837 55 77.49633 84.57396    bc  1.626823
#Comparing the levels of A for each level of B
ggplot(ABC.means.cld, aes(x = A, 
                         y = emmean,
                         fill = B, 
                         label = .group)) + 
  geom_col(position = "dodge") +
  geom_errorbar(aes(ymin = emmean - SE1, 
                    ymax = emmean + SE1),
                width = 0.2, 
                size = 0.7,
                position = position_dodge(0.8)) +
  geom_text(vjust=-2.5, 
            position = position_dodge(0.9),
            size = 2.5) +
  labs(x= "A", y= "Mean Yield") +
  lims(y=c(0,100))+
  facet_wrap(. ~ C) +
  theme_classic() +
  theme(legend.position = "bottom") +
   theme(legend.title=element_blank())