Graeco-Latin Square Design

STA6246

Author

Dr. Cohen

Introduction

  • The Graeco-Latin Square Design (GLSD) allows us to eliminate three nuisance factors (blocking).

  • Each factor has p levels

  • We can control the effect of these blocking factors by grouping the experimental units into blocks.

  • This is an extension of the Latin Square Design.

Rocket Propellant Experiment

  • An experimenter wanted to study the effects of 5 different formulations of a rocket propellant on the observed burning rate.

  • Each formulation is mixed from a batch of raw material.

  • Each formulation is prepared by one operator for each batch.

  • Each formulation can be assembly tested using 5 different tests.

  • \(p=5\) levels for each factor

\(p^2= 25\) runs or observations.

Data description

library(tidyverse)
library(gtsummary)
# enter ER
BurnRate=c(24,17,18,26,22,20,24,38,31,30,19,30,26,26,20,24,27,27,23,29,24,36,21,22,31)
Bt=gl(n = 5, k = 1,length = 25, labels = c(1:5))
Op=gl(n = 5, k = 5,length = 25, labels = c(1:5))
Trts=c("A","B","C","D","E",
       "B","C","D","E","A",
       "C","D","E","A","B",
       "D","E","A","B","C",
       "E","A","B","C","D")

# ADDITIONAL FACTOR - THAT IS THE TEST ASSEMBLIES
Tests=c("Aa","Ba","Ca","Da","Ea",
        "Da","Ea","Aa","Ba","Ca",
        "Ba","Ca","Da","Ea","Aa",
        "Ea","Aa","Ba","Ca","Da",
        "Ca","Da","Ea","Aa","Ba")

rocket=tibble(BurnRate,Bt,Op,Trts,Tests)

head(rocket)
# A tibble: 6 × 5
  BurnRate Bt    Op    Trts  Tests
     <dbl> <fct> <fct> <chr> <chr>
1       24 1     1     A     Aa   
2       17 2     1     B     Ba   
3       18 3     1     C     Ca   
4       26 4     1     D     Da   
5       22 5     1     E     Ea   
6       20 1     2     B     Da   
rocket %>%
  ggplot(aes(x=Trts,y=BurnRate,group=Trts)) +
  geom_point()+
  stat_summary(aes(group=1),
               fun = mean, 
               geom = "line", 
               col="red",
               linewidth=2)+
  stat_summary(aes(label=..y..,group=1),
               fun = mean, 
               geom = "text", 
               col="blue",
               size=10)+
  theme_bw()+
  labs(x="Formualtions",
       y="Burn Rate")

rocket %>% 
  select(BurnRate,Trts) %>%
  tbl_summary(by=Trts,
              statistic = list(all_continuous() ~ "{mean} ({sd})"))
Characteristic A, N = 51 B, N = 51 C, N = 51 D, N = 51 E, N = 51
BurnRate 28.6 (4.7) 20.2 (2.2) 22.4 (4.4) 29.8 (5.4) 26.0 (3.4)
1 Mean (SD)

  • The mean burning rate seems to vary acrross the formulations.

  • The variability seems to vary a bit across the formulations

  • The Formulation D provides the highest mean burning rate and B the lowest.

ANOVA in R

Fitting
library(sjPlot)

# ANOVA one factor and one blocking factor
results = aov(BurnRate ~ Trts + Bt + Op + Tests ,data = rocket)
summary(results)
            Df Sum Sq Mean Sq F value Pr(>F)   
Trts         4  330.0    82.5   7.933 0.0069 **
Bt           4   68.0    17.0   1.635 0.2566   
Op           4  150.0    37.5   3.606 0.0579 . 
Tests        4   44.8    11.2   1.077 0.4284   
Residuals    8   83.2    10.4                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
tab_model(results)
  BurnRate
Predictors p
Trts 0.007
Bt 0.257
Op 0.058
Tests 0.428
Residuals
Observations 25
R2 / R2 adjusted 0.877 / 0.631 
# Individual CIs 
confint(results)
                 2.5 %    97.5 %
(Intercept)  16.067590 28.332410
TrtsB       -13.103344 -3.696656
TrtsC       -10.903344 -1.496656
TrtsD        -3.503344  5.903344
TrtsE        -7.303344  2.103344
Bt2          -0.103344  9.303344
Bt3          -0.903344  8.503344
Bt4          -1.303344  8.103344
Bt5          -0.503344  8.903344
Op2           2.496656 11.903344
Op3          -1.903344  7.503344
Op4          -0.103344  9.303344
Op5           0.696656 10.103344
TestsBa      -5.903344  3.503344
TestsCa      -5.903344  3.503344
TestsDa      -3.503344  5.903344
TestsEa      -7.503344  1.903344

The ANOVA was rejected (p-value \(= 0.007\)). Therefore, the data is compatible with the means difference in the burn rate across the 5 formulations. Next, we will diagnose the model.

Model Adequacy
par(mfrow=c(2,2))
plot(results)

We could consider a transformation because the Scale-Location plots is showing a trend.

Transformation

Since \(\lambda\) was estimated to be \(0\) we need a \(\log\) transformation.

library(car)
l=powerTransform(results)
l$roundlam
Y1 
 0 
## do the transofrmation
Tburn =log(BurnRate)
## Refit ANOVA
rano=aov(Tburn ~ Bt + Op + Trts + Tests, data=rocket)
## Summary
summary(rano)
            Df Sum Sq Mean Sq F value  Pr(>F)   
Bt           4 0.0933 0.02333   1.870 0.20912   
Op           4 0.2316 0.05789   4.641 0.03122 * 
Trts         4 0.5301 0.13252  10.624 0.00275 **
Tests        4 0.0588 0.01469   1.178 0.38956   
Residuals    8 0.0998 0.01247                   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
par(mfrow=c(2,2))
plot(rano)

Multiple Comparison

Considering 5% level of significance.

#Tukey Test 
THSD=TukeyHSD(results,which = "Trts")
plot(THSD,las=1)

THSD
  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = BurnRate ~ Trts + Bt + Op + Tests, data = rocket)

$Trts
    diff        lwr       upr     p adj
B-A -8.4 -15.446335 -1.353665 0.0205953
C-A -6.2 -13.246335  0.846335 0.0881536
D-A  1.2  -5.846335  8.246335 0.9731074
E-A -2.6  -9.646335  4.446335 0.7124698
C-B  2.2  -4.846335  9.246335 0.8126378
D-B  9.6   2.553665 16.646335 0.0097262
E-B  5.8  -1.246335 12.846335 0.1152305
D-C  7.4   0.353665 14.446335 0.0395290
E-C  3.6  -3.446335 10.646335 0.4510964
E-D -3.8 -10.846335  3.246335 0.4045667
# Fisher LSD Test
library(agricolae)
L=LSD.test(results,"Trts",console=TRUE,group = F)

Study: results ~ "Trts"

LSD t Test for BurnRate 

Mean Square Error:  10.4 

Trts,  means and individual ( 95 %) CI

  BurnRate      std r      LCL      UCL Min Max
A     28.6 4.669047 5 25.27423 31.92577  24  36
B     20.2 2.167948 5 16.87423 23.52577  17  23
C     22.4 4.393177 5 19.07423 25.72577  18  29
D     29.8 5.403702 5 26.47423 33.12577  24  38
E     26.0 3.391165 5 22.67423 29.32577  22  31

Alpha: 0.05 ; DF Error: 8
Critical Value of t: 2.306004 

Comparison between treatments means

      difference pvalue signif.        LCL       UCL
A - B        8.4 0.0034      **   3.696656 13.103344
A - C        6.2 0.0161       *   1.496656 10.903344
A - D       -1.2 0.5725          -5.903344  3.503344
A - E        2.6 0.2382          -2.103344  7.303344
B - C       -2.2 0.3122          -6.903344  2.503344
B - D       -9.6 0.0015      ** -14.303344 -4.896656
B - E       -5.8 0.0217       * -10.503344 -1.096656
C - D       -7.4 0.0067      ** -12.103344 -2.696656
C - E       -3.6 0.1156          -8.303344  1.103344
D - E        3.8 0.0995       .  -0.903344  8.503344

Conclusion:

  • Tukey’s Test indicates:

    • Significant differences between:

      • B and D;

      • A and B; and

      • C and D.

  • Fisher LSD Test indicates:

    • Significant differences between:

      • B and D;

      • A and B;

      • C and D;

      • B and E;

      • A and C.

Formulation D led to the highest mean burning rate and it is statistical significant different from formations B and C.

Formulation D, A, and E resulted in no significant differences.