1 Latin-Square

Author:

- Author 1 Long Vu Tran

- Author 2 A K M Sarower Kabir

- Author 3 Manoj Giri

The Latin square design is used to eliminate two known and controllable nuisance sources of variability; that is, it systematically allows blocking in two directions.

The number of levels of both blocking variables needs to be the same and must equal the number of treatment levels.

Linear effect equation:

\[ Y_{ijk}=\mu+\tau_{i}+\beta_{j}+\alpha_{k}+\epsilon_{ijk} \]

Where,

\(\mu\) = Grand Mean

\(\tau_{i}\) = Treatment effect

\(\beta_{j}\) = Block-1 effect

\(\alpha_{k}\) = Block-2 effect

\(\epsilon_{ijk}\) = Random error

i = number of treatments

j = number of block-1

k = number of block-2

1.1 Assumptions

The errors \(\epsilon_{ijk}\) are normally and independently distributed with mean zero and constant variance \(\sigma^{2}\).
There is no interaction between the two blocks and treatment.

1.2 Process

1.2.1 Sample Size Determination

For a p x p Latin square, the number of observations for each replication is \(p^{2}\).

When small Latin squares are used, it is frequently desirable to replicate them to increase the error degrees of freedom. To determine the number of replications, estimate the power of the test for each number of replications then compare it with the expected power. There are two cases to estimate the power of the test:

Pre-sample available with known variation

#power.anova.test(groups= 'number of group',n= 'sample size',between.var= 'variation between group',within.var= 'variation within group',sig.level= alpha)

Without pre-samples

#d is the range of means divided by sigma, k is the number of populations
#without pre-samples, d is specified as the difference in means as a proportion of sigma to which the power correspondes
#The effect is given by d*sqrt(1/(2*k)) 
library(pwr)
#pwr.anova.test(k='number of group',n='sample size',f=d*sqrt(1/(2*k)),sig.level=alpha)
#plot(pwr.anova.test(k='number of group',n=sample size',d*sqrt(1/(2*k)),sig.level=alpha)

1.2.2 Design Layout

Use the function OutDesign below to export the “Latin square layout.csv” file in the working directory folder.

library(agricolae)
OutDesign<-function(treatment,block1,block2){
  outdesign<-design.lsd(treatment, )$book
  row <- seq(1,length(treatment))
  outdesign <- merge(outdesign, data.frame(block1,row), by = "row")
  col <- seq(1,length(treatment))
  outdesign <- merge(outdesign,data.frame(block2, col), by = 'col')
  outdesign <- outdesign[,3:6]
  write.csv(outdesign,"Latin square layout.csv")}

1.2.3 Collect Data

Adding the response column then conduct the experiment to full fill that column.

Use the InDesign function below to input back the file into R.

InDesign<-function(path){
  indesign<-read.csv(path)
  indesign<-indesign[,3:6]
  indesign$Response <- as.numeric(indesign$Response)
  indesign
  }

1.2.4 Preliminary Plots

A histogram of response data can be used to have an initial look on the distribution.

#hist("Response, main="Title of the Histogram", Xlab="X label", ylab="Y label")

A Box-plot also can be used to review if there are outliers in response.

#boxplot('Response')

1.2.5 Statistical Test

Null hypotheses: \(H_{o}:\mu_{ 1}=\mu_{ 2}= ... =\mu_{ i}\)

Alternative Hypotheses: \(H_{a}:\) At least one \(\mu\) is different

Create the ANOVA table:

#model <- aov(Response~treatment+Block1+Block2)
#summary(model)

1.2.6 Residuals

Create plots of residuals to check the assumption:

#plot(model)

Normal probability plot: If the underlying error distribution is normal, this plot will resemble a straight line.
Plot of Residuals Versus Fitted Values: the residuals should be structureless and unrelated to the predicted response.
Plots of Residuals Versus treatment levels: Variation of residuals should be equal for all levels of treatment.

When the assumptions are violated, Variance Stabilizing Transformation (VST) using Box-Cox or Non-Parametric analysis should be used.

1.2.7 Inference

If the assumptions are met, then compare the P-value of treatment in the ANOVA table with \(\alpha\) to have the conclusion on the effect of treatment.

1.2.8 Multiple Comparisons

After rejecting the \(H_{o}\) hypothesis, the two most commonly used techniques for making multiple comparisons between group means are LSD (Least Significant Difference) and HSD (Honestly Significant Difference) or also called Tukey’s HSD Test.

LSD test (Least Significant Difference test):
1)Performs all possible pairwise comparisons between group means.
2)Uses the pooled standard error from the ANOVA to calculate a t-statistic and p-value for each pair of means.

library(agricolae)
##LSD.test(y, trt, Dferror, MSerror, alpha)     
#y= model(aov or lm) or answer of the experimental unit     
#trt= Constant( only y=model) or vector treatment applied to each experimental unit.    
#DFerror= Degrees of freedom of the experimental error.    
#MSerror= Means square error of the experimental.    
#alpha= Level of risk for the test

If we select \(\alpha\) = 0.05, the probability of reaching the correct decision on any single comparison is 0.95. However, the probability of reaching the correct conclusion on all comparisons is considerably less than 0.95, so the type I error is inflated in LSD test. Whereas TukeyHSD will control the type 1 error rate by adjusting the p-value.

Tukey’s HSD Test:
1)Uses the Studentized range statistic to make all pairwise comparisons.
2)Controls family-wise error rate for multiple testing.

library(car)
##TukeyHSD(x)    
#x=A fitted model object, usually an aov fit.     
##plot(TukeyHSD())

1.3 Example

Automobile Emission Experiment.

Four cars and four drivers are employed in a study for possible differences in effect between four gasoline additives(A, B, C, D) on the emission of the car. Even though cars can be identical models, slight systematic differences are likely to occur in their performance, and even though each driver may do his best to drive the car in the manner required by the test, slight systematic differences can occur from driver to driver. It would be desirable to eliminate both the car-to-car and driver-to-driver differences. With the \(\alpha\) = 0.05, conduct the experiment and then analyze the effect of gasoline additives on emission.

1.3.1 Sample size determination:

Since the gasoline additives have 4 levels (A, B, C, D), the sample size of Latin square design is 4 x 4 = 16 for one replication.

1.3.2 Create the design layout:

GasAdditives <- c('A','B','C','D')
Car <- c('Car1', 'Car2', 'Car3', 'Car4')
Driver <- c('Driver1', 'Driver2', 'Driver3', 'Driver4')
OutDesign(GasAdditives, Car, Driver)

After collecting data, use the InDesign function to import into R

indesign <- InDesign('https://raw.githubusercontent.com/tlvu-mecha-data/DOEcookbookLatinsquare/main/dataLatinsquarelayout.txt')

1.3.3 Preliminary plots of the response:

Histogram:

hist(indesign$Response, main="Histogram of Emission", xlab="Emission", col ='blue')

Box-plot:

boxplot(Response~treatment, data = indesign, main = 'Box-plot of Emission', ylab = 'Emission', xlab = 'Treatment', col = 'green')

The variance on treatment A are higher than B, C and D.
Treatment D has the lowest median of emission.

1.3.4 Statistical test:

Hypothesis:

Null hypotheses: \(H_{o}:\mu_{ 1}=\mu_{ 2}= \mu_{3} =\mu_{ 4}\)

Alternative Hypotheses: \(H_{a}:\) At least one \(\mu\) is different

Create the ANOVA table:

model <- aov(indesign$Response~indesign$treatment+indesign$block1+indesign$block2)
summary(model)

##                    Df Sum Sq Mean Sq F value   Pr(>F)    
## indesign$treatment  3 182.75   60.92  27.074 0.000693 ***
## indesign$block1     3  38.25   12.75   5.667 0.034812 *  
## indesign$block2     3   1.25    0.42   0.185 0.902711    
## Residuals           6  13.50    2.25                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

1.3.5 Residuals

plot(model)

Normal probability plot: This plot resemble a straight line indicate the residuals follow normal distribution.
Plot of Residuals Versus Fitted Values: the residuals are structureless and unrelated to the predicted response.
Plots of Residuals Versus treatment levels: Variation of residuals are equal for all levels of treatment.

1.3.6 Inference:

The assumptions are met, P-value = 0.000693 which is smaller than \(\alpha\) = 0.05.

=> Reject the Null hypothesis

1.3.7 Multiple Comparisons

LSD test (Least Significant Difference test):

library(agricolae)
indesign$treatment <- as.factor(indesign$treatment)
lsd <- LSD.test(y = indesign$Response, trt = indesign$treatment, DFerror = 6, MSerror = 2.25, alpha = 0.05)
lsd

## $statistics
##   MSerror Df   Mean       CV  t.value      LSD
##      2.25  6 22.125 6.779661 2.446912 2.595342
## 
## $parameters
##         test p.ajusted             name.t ntr alpha
##   Fisher-LSD      none indesign$treatment   4  0.05
## 
## $means
##   indesign$Response      std r   se      LCL      UCL Min Max   Q25  Q50   Q75
## A             20.75 2.629956 4 0.75 18.91482 22.58518  18  23 18.75 21.0 23.00
## B             25.75 1.707825 4 0.75 23.91482 27.58518  24  28 24.75 25.5 26.50
## C             24.75 2.217356 4 0.75 22.91482 26.58518  22  27 23.50 25.0 26.25
## D             17.25 1.707825 4 0.75 15.41482 19.08518  15  19 16.50 17.5 18.25
## 
## $comparison
## NULL
## 
## $groups
##   indesign$Response groups
## B             25.75      a
## C             24.75      a
## A             20.75      b
## D             17.25      c
## 
## attr(,"class")
## [1] "group"

plot(lsd, main = 'LSD chart for treatment', xlab = 'Gasoline additive', ylab = 'Emission')

Comments:

From the LSD test mean of treatment B and C are equal. Mean of treatment A is different from B, C and D. Mean of treatment D is different from B, C and A.

Tukey’s HSD Test:

Tukey <- TukeyHSD(model, conf.level = 0.95)
Tukey

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = indesign$Response ~ indesign$treatment + indesign$block1 + indesign$block2)
## 
## $`indesign$treatment`
##     diff         lwr        upr     p adj
## B-A  5.0   1.3283006  8.6716994 0.0129205
## C-A  4.0   0.3283006  7.6716994 0.0351677
## D-A -3.5  -7.1716994  0.1716994 0.0603481
## C-B -1.0  -4.6716994  2.6716994 0.7845973
## D-B -8.5 -12.1716994 -4.8283006 0.0008344
## D-C -7.5 -11.1716994 -3.8283006 0.0016471
## 
## $`indesign$block1`
##            diff         lwr       upr     p adj
## Car2-Car1  3.75  0.07830061 7.4216994 0.0459307
## Car3-Car1  0.75 -2.92169939 4.4216994 0.8908002
## Car4-Car1  3.00 -0.67169939 6.6716994 0.1056937
## Car3-Car2 -3.00 -6.67169939 0.6716994 0.1056937
## Car4-Car2 -0.75 -4.42169939 2.9216994 0.8908002
## Car4-Car3  2.25 -1.42169939 5.9216994 0.2472640
## 
## $`indesign$block2`
##                  diff       lwr      upr     p adj
## Driver2-Driver1 -0.25 -3.921699 3.421699 0.9948966
## Driver3-Driver1  0.25 -3.421699 3.921699 0.9948966
## Driver4-Driver1 -0.50 -4.171699 3.171699 0.9626993
## Driver3-Driver2  0.50 -3.171699 4.171699 0.9626993
## Driver4-Driver2 -0.25 -3.921699 3.421699 0.9948966
## Driver4-Driver3 -0.75 -4.421699 2.921699 0.8908002

plot(Tukey, col = 'blue')