Data

The effect of five ingredients (A, B, C, D, E) on the reaction time of a chemical process is studied. Each batch of material allows only five runs, and only five runs can be made in one day, so the experiment is run as a Latin square with Batch as the rows and Day as the columns. The data are entered in the same wide format as the table (each cell is “ingredient = reaction time”), converted to a tidy format using pivot_longer, and then split into the ingredient letter and the reaction time using separate.

  Batch Day Ingredient Time
1     1   1          A    8
2     1   2          B    7
3     1   3          D    1
4     1   4          C    7
5     1   5          E    3
6     2   1          C   11
7     2   2          E    2
8     2   3          A    7
  Ingredient Mean Variance
1          A  8.4      1.3
2          B  5.6      4.3
3          C  8.8      2.7
4          D  3.4      4.3
5          E  3.2      3.7
  Batch Mean Variance
1     1  5.2      9.2
2     2  6.2     13.7
3     3  5.8     13.7
4     4  7.2      3.2
5     5  5.0      8.0
  Day Mean Variance
1   1  6.6      8.8
2   2  5.6     11.3
3   3  5.4     12.3
4   4  5.0      8.5
5   5  6.8      7.7

Question 1

Yes, this is a valid Latin square. A p × p Latin square arranges p treatments in p rows and p columns so that each treatment appears exactly once in every row and exactly once in every column. Here p = 5: there are 5 ingredients (A–E), 5 batches (rows) and 5 days (columns), and every batch–day cell contains exactly one run, for 25 runs in total.

The tables below count how many times each ingredient appears in each batch and in each day.

     Ingredient
Batch A B C D E
    1 1 1 1 1 1
    2 1 1 1 1 1
    3 1 1 1 1 1
    4 1 1 1 1 1
    5 1 1 1 1 1
   Ingredient
Day A B C D E
  1 1 1 1 1 1
  2 1 1 1 1 1
  3 1 1 1 1 1
  4 1 1 1 1 1
  5 1 1 1 1 1
[1] TRUE

Every count is exactly 1, so each ingredient is used once in every batch and once on every day. Because of this balance, each ingredient mean is averaged over all five batches and all five days, so batch-to-batch and day-to-day differences cancel out of the comparison between ingredients. The square is not in standard form (the first row is A, B, D, C, E and the first column is A, C, B, D, E, not alphabetical), which is fine; it is simply one of the many randomized 5 × 5 Latin squares.

Question 2

The linear effects model for the Latin square design is

\[y_{ijk} = \mu + \alpha_i + \tau_j + \beta_k + \varepsilon_{ijk}, \qquad i = 1, 2, \ldots, 5, \quad j = 1, 2, \ldots, 5, \quad k = 1, 2, \ldots, 5\]

where \(y_{ijk}\) is the reaction time in batch \(i\) with ingredient \(j\) on day \(k\), \(\mu\) is the overall mean, \(\alpha_i\) is the effect of the \(i\)-th batch (row), \(\tau_j\) is the effect of the \(j\)-th ingredient (treatment), \(\beta_k\) is the effect of the \(k\)-th day (column), and \(\varepsilon_{ijk} \sim NID(0, \sigma^2)\) is the random error. The effects satisfy \(\sum_{i} \alpha_i = \sum_{j} \tau_j = \sum_{k} \beta_k = 0\).

Only 25 of the 125 possible \((i, j, k)\) combinations are observed, because once the batch \(i\) and the day \(k\) are chosen, the ingredient \(j\) is fixed by the square. The model is additive, so it assumes there is no interaction between batches, ingredients and days.

The hypotheses being tested for the factor of interest (Ingredient) are:

Question 3

All three factors are declared as factors (as.factor()), and the model is fit with aov().

            Df Sum Sq Mean Sq F value   Pr(>F)    
Ingredient   4 141.44   35.36  11.309 0.000488 ***
Batch        4  15.44    3.86   1.235 0.347618    
Day          4  12.24    3.06   0.979 0.455014    
Residuals   12  37.52    3.13                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Ingredient: F₀ = 11.31 with a p-value of 0.000488. Since the p-value (0.000488) < α = 0.05, we reject H₀ and conclude that the ingredients have a significant effect on the reaction time of the chemical process.

Batch and Day: F₀ = 1.235 (p-value = 0.3476) for batch and F₀ = 0.979 (p-value = 0.4550) for day. Neither nuisance factor is significant at α = 0.05; their mean squares (3.86 and 3.06) are about the same size as the error mean square (3.13), so the batch-to-batch and day-to-day variability turned out to be small in this experiment.

To see which ingredients differ, Tukey’s test is applied to the ingredient means.

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

Fit: aov(formula = Time ~ Ingredient + Batch + Day, data = dat)

$Ingredient
    diff        lwr        upr     p adj
B-A -2.8 -6.3646078  0.7646078 0.1539433
C-A  0.4 -3.1646078  3.9646078 0.9960012
D-A -5.0 -8.5646078 -1.4353922 0.0055862
E-A -5.2 -8.7646078 -1.6353922 0.0041431
C-B  3.2 -0.3646078  6.7646078 0.0864353
D-B -2.2 -5.7646078  1.3646078 0.3365811
E-B -2.4 -5.9646078  1.1646078 0.2631551
D-C -5.4 -8.9646078 -1.8353922 0.0030822
E-C -5.6 -9.1646078 -2.0353922 0.0023007
E-D -0.2 -3.7646078  3.3646078 0.9997349

Four pairs are significantly different at α = 0.05: D–A, E–A, D–C and E–C (all p-values < 0.01). Ingredient B is not significantly different from any other ingredient. So ingredients A and C (means 8.4 and 8.8) give the longest reaction times, ingredients D and E (means 3.4 and 3.2) give the shortest, and B (mean 5.6) lies in between.

The plot shows the same picture: A and C sit high, D and E sit low, and B is in the middle.

The model assumptions are checked with the residuals.


    Shapiro-Wilk normality test

data:  residuals(model_ls)
W = 0.96606, p-value = 0.5476

The residuals versus fitted values show no clear pattern or funnel shape, and the normal Q-Q plot is close to a straight line (Shapiro–Wilk p-value = 0.5476), so the normality and constant-variance assumptions are reasonable and the ANOVA results can be trusted.

Conclusion: The ingredient has a significant effect on reaction time (p = 0.000488). Ingredients D and E give the shortest reaction times, and A and C the longest. If a short reaction time is desired, ingredient D or E should be used.

R Code

library(dplyr)
library(tidyr)
library(ggplot2)

# Data in tidy format
wide <- data.frame(
  Batch = 1:5,
  Day_1 = c("A=8",  "C=11", "B=4",  "D=6",  "E=4"),
  Day_2 = c("B=7",  "E=2",  "A=9",  "C=8",  "D=2"),
  Day_3 = c("D=1",  "A=7",  "C=10", "E=6",  "B=3"),
  Day_4 = c("C=7",  "D=3",  "E=1",  "B=6",  "A=8"),
  Day_5 = c("E=3",  "B=8",  "D=5",  "A=10", "C=8")
)

dat <- wide %>%
  pivot_longer(cols = -Batch, names_to = "Day", names_prefix = "Day_", values_to = "Cell") %>%
  separate(Cell, into = c("Ingredient", "Time"), sep = "=", convert = TRUE) %>%
  mutate(Batch = as.factor(Batch),
         Day = as.factor(Day),
         Ingredient = as.factor(Ingredient))

head(as.data.frame(dat), 8)

dat %>%
  group_by(Ingredient) %>%
  summarise(Mean = round(mean(Time), 4), Variance = round(var(Time), 4)) %>%
  as.data.frame()

dat %>%
  group_by(Batch) %>%
  summarise(Mean = round(mean(Time), 4), Variance = round(var(Time), 4)) %>%
  as.data.frame()

dat %>%
  group_by(Day) %>%
  summarise(Mean = round(mean(Time), 4), Variance = round(var(Time), 4)) %>%
  as.data.frame()

# Question 1: check that this is a valid Latin square
table(Batch = dat$Batch, Ingredient = dat$Ingredient)
table(Day = dat$Day, Ingredient = dat$Ingredient)
all(table(dat$Batch, dat$Ingredient) == 1) & all(table(dat$Day, dat$Ingredient) == 1)

# Question 3: Latin square ANOVA with aov()
model_ls <- aov(Time ~ Ingredient + Batch + Day, data = dat)
summary(model_ls)

# Which ingredients differ
TukeyHSD(model_ls, "Ingredient", conf.level = 0.95)

ggplot(dat, aes(x = Ingredient, y = Time)) +
  geom_point(size = 2, color = "grey50") +
  stat_summary(fun = mean, geom = "point", size = 4, color = "red") +
  stat_summary(fun = mean, geom = "line", aes(group = 1), color = "red") +
  labs(x = "Ingredient", y = "Reaction Time",
       title = "Reaction Time by Ingredient (red = ingredient mean)")

# Residual checks
par(mfrow = c(1, 2))
plot(model_ls, which = 1)
plot(model_ls, which = 2)
par(mfrow = c(1, 1))

shapiro.test(residuals(model_ls))