Latin Square Design

STA6246 - DOE

Author

Dr. Cohen

Introduction

  • The Latin Square Design (LSD) allows us to eliminate two nuisance factors (blocking).
  • Each factor has p levels
  • We can control the effect of these blocking factors by grouping the experimental units into blocks.
  • We consider a square wit p rows and p columns.

We can geenrate latin square in R as follows:

library(magic)
# generate r by r latin square design
rlatin(5)
     [,1] [,2] [,3] [,4] [,5]
[1,]    4    2    1    5    3
[2,]    5    1    3    2    4
[3,]    1    5    4    3    2
[4,]    3    4    2    1    5
[5,]    2    3    5    4    1

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.

  • \(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")

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

rocket %>%
  ggplot(aes(x=Trts,y=BurnRate)) +
  geom_point()+
  theme_bw()+
  labs(x="Formualtions",
       y="Burn Rate")

rocket %>%
  tbl_summary(by=Trts)
Characteristic A, N = 51 B, N = 51 C, N = 51 D, N = 51 E, N = 51
BurnRate 27.0 (26.0, 30.0) 20.0 (20.0, 21.0) 22.0 (19.0, 24.0) 30.0 (26.0, 31.0) 26.0 (24.0, 27.0)
Bt
    1 1 (20%) 1 (20%) 1 (20%) 1 (20%) 1 (20%)
    2 1 (20%) 1 (20%) 1 (20%) 1 (20%) 1 (20%)
    3 1 (20%) 1 (20%) 1 (20%) 1 (20%) 1 (20%)
    4 1 (20%) 1 (20%) 1 (20%) 1 (20%) 1 (20%)
    5 1 (20%) 1 (20%) 1 (20%) 1 (20%) 1 (20%)
Op
    1 1 (20%) 1 (20%) 1 (20%) 1 (20%) 1 (20%)
    2 1 (20%) 1 (20%) 1 (20%) 1 (20%) 1 (20%)
    3 1 (20%) 1 (20%) 1 (20%) 1 (20%) 1 (20%)
    4 1 (20%) 1 (20%) 1 (20%) 1 (20%) 1 (20%)
    5 1 (20%) 1 (20%) 1 (20%) 1 (20%) 1 (20%)
1 Median (IQR); n (%)

  • The formulations seem to impact the burnig rate

  • The variability seems to be similar across the formulations

ANOVA in R

Fitting
library(sjPlot)

# ANOVA one factor and one blocking factor
results = aov(BurnRate ~ Trts + Bt + Op - 1 ,data = rocket)
summary(results)
          Df Sum Sq Mean Sq F value   Pr(>F)    
Trts       5  16459    3292 308.606 3.07e-12 ***
Bt         4     68      17   1.594   0.2391    
Op         4    150      38   3.516   0.0404 *  
Residuals 12    128      11                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
tab_model(results)
  BurnRate
Predictors p
Trts <0.001
Bt 0.239
Op 0.040
Residuals
Observations 25
R2 / R2 adjusted 0.992 / 0.984 
# Individual CIs - need to add the CI for the reference level
confint(results)
            2.5 %    97.5 %
TrtsA 16.26859896 26.531401
TrtsB  7.86859896 18.131401
TrtsC 10.06859896 20.331401
TrtsD 17.46859896 27.731401
TrtsE 13.66859896 23.931401
Bt2    0.09946357  9.100536
Bt3   -0.70053643  8.300536
Bt4   -1.10053643  7.900536
Bt5   -0.30053643  8.700536
Op2    2.69946357 11.700536
Op3   -1.70053643  7.300536
Op4    0.09946357  9.100536
Op5    0.89946357  9.900536

The ANOVA was rejected (p-value \(= 0.003\)). Therefore, the data is compatible with the means difference in the burn rate across the 5 fomulations. 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
library(car)
l=powerTransform(results)
l$roundlam
Y1 
 0 
Tburn =log(BurnRate)
rano=aov(Tburn ~ Bt + Op + Trts)
summary(rano)
            Df Sum Sq Mean Sq F value   Pr(>F)    
Bt           4 0.0933 0.02333   1.766 0.200512    
Op           4 0.2316 0.05789   4.382 0.020557 *  
Trts         4 0.5301 0.13252  10.030 0.000837 ***
Residuals   12 0.1586 0.01321                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
boxplot(Tburn~Trts)

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)

# 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.66667 

Trts,  means and individual ( 95 %) CI

  BurnRate      std r      LCL      UCL Min Max
A     28.6 4.669047 5 25.41764 31.78236  24  36
B     20.2 2.167948 5 17.01764 23.38236  17  23
C     22.4 4.393177 5 19.21764 25.58236  18  29
D     29.8 5.403702 5 26.61764 32.98236  24  38
E     26.0 3.391165 5 22.81764 29.18236  22  31

Alpha: 0.05 ; DF Error: 12
Critical Value of t: 2.178813 

Comparison between treatments means

      difference pvalue signif.         LCL        UCL
A - B        8.4 0.0016      **   3.8994636 12.9005364
A - C        6.2 0.0110       *   1.6994636 10.7005364
A - D       -1.2 0.5720          -5.7005364  3.3005364
A - E        2.6 0.2321          -1.9005364  7.1005364
B - C       -2.2 0.3078          -6.7005364  2.3005364
B - D       -9.6 0.0006     *** -14.1005364 -5.0994636
B - E       -5.8 0.0158       * -10.3005364 -1.2994636
C - D       -7.4 0.0038      ** -11.9005364 -2.8994636
C - E       -3.6 0.1069          -8.1005364  0.9005364
D - E        3.8 0.0907       .  -0.7005364  8.3005364

Conclusion:

Formation D led to the higher mean burning rate and it is statistical significant different from formation C (p=0.025; 95% CI 0.82 to 14) and B (p=0.004; 95% CI 3.016 to 16.18).