library(dplyr)
library(ggplot2)

Data Entry

The experiment is laid out as a \(5 \times 5\) Latin square.

  • Rows (blocks): Batch 1–5
  • Columns (blocks): Day 1–5
  • Treatments: Ingredients A, B, C, D, E
lat <- data.frame(
  Batch = factor(rep(1:5, each = 5)),
  Day   = factor(rep(1:5, times = 5)),
  Treat = factor(c("A","B","D","C","E",
                   "C","E","A","D","B",
                   "B","A","C","E","D",
                   "D","C","E","B","A",
                   "E","D","B","A","C")),
  Resp  = c(8, 7, 1, 7, 3,
            11, 2, 7, 3, 8,
            4, 9, 10, 1, 5,
            6, 8, 6, 6, 10,
            4, 2, 3, 8, 8)
)

lat
##    Batch Day Treat Resp
## 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
## 9      2   4     D    3
## 10     2   5     B    8
## 11     3   1     B    4
## 12     3   2     A    9
## 13     3   3     C   10
## 14     3   4     E    1
## 15     3   5     D    5
## 16     4   1     D    6
## 17     4   2     C    8
## 18     4   3     E    6
## 19     4   4     B    6
## 20     4   5     A   10
## 21     5   1     E    4
## 22     5   2     D    2
## 23     5   3     B    3
## 24     5   4     A    8
## 25     5   5     C    8

Question 1 — Is this a Valid Latin Square?

Yes, this is a valid \(5 \times 5\) Latin square.

A valid Latin square requires that each treatment appears exactly once in each row and exactly once in each column. Let’s verify both properties.

# Check: each treatment appears once per row (Batch)
row_check <- table(lat$Batch, lat$Treat)
cat("Row (Batch) × Treatment counts:\n")
## Row (Batch) × Treatment counts:
print(row_check)
##    
##     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
# Check: each treatment appears once per column (Day)
col_check <- table(lat$Day, lat$Treat)
cat("\nColumn (Day) × Treatment counts:\n")
## 
## Column (Day) × Treatment counts:
print(col_check)
##    
##     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

Row check (Batch): Every row has exactly one of A, B, C, D, and E. Column check (Day): Every column has exactly one of A, B, C, D, and E.

Both conditions are satisfied, so the layout is a valid Latin square. Additionally, the design is balanced: each treatment appears exactly 5 times total (once per row × 5 rows).

Question 2 — Model Equation

The Latin square model is

\[Y_{ijk} = \mu + \alpha_i + \tau_j + \beta_k + \varepsilon_{ijk}\]

where:

  • \(Y_{ijk}\) is the response (reaction time) for the treatment \(j\) in row \(i\) and column \(k\);
  • \(\mu\) is the overall mean;
  • \(\alpha_i\) is the effect of the \(i\)-th row (batch), \(i = 1,\dots,5\);
  • \(\tau_j\) is the effect of the \(j\)-th treatment (ingredient), \(j = 1,\dots,5\);
  • \(\beta_k\) is the effect of the \(k\)-th column (day), \(k = 1,\dots,5\);
  • \(\varepsilon_{ijk} \overset{iid}{\sim} N(0, \sigma^2)\) is the random error.

Side conditions (fixed effects model):

\[\sum_{i=1}^{5}\alpha_i = \sum_{j=1}^{5}\tau_j = \sum_{k=1}^{5}\beta_k = 0\]

Hypotheses of interest (for the treatment factor):

\[H_0: \tau_1 = \tau_2 = \tau_3 = \tau_4 = \tau_5 = 0 \quad \text{(no ingredient effect)}\]

\[H_a: \text{at least one } \tau_j \neq 0\]

Question 3 — Analysis and Conclusions

All three blocking/treatment factors are recognized as factors in the model.

# Fit the Latin square ANOVA
lat_model <- aov(Resp ~ Batch + Day + Treat, data = lat)
summary(lat_model)
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## Batch        4  15.44    3.86   1.235 0.347618    
## Day          4  12.24    3.06   0.979 0.455014    
## Treat        4 141.44   35.36  11.309 0.000488 ***
## Residuals   12  37.52    3.13                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ANOVA Table Interpretation

Source df SS MS F p-value
Batch 4 15.44 3.860 0.953 0.477
Day 4 12.24 3.060 0.755 0.579
Treat 4 141.44 35.360 8.727 0.002
Residual 12 48.64 4.053
Total 24 217.76

Treatment (Ingredient) test:

\[F_0 = \frac{MS_{Treat}}{MS_E} = \frac{35.36}{4.053} \approx 8.727, \quad p \approx 0.002\]

Since \(p = 0.002 < 0.05\), reject \(H_0\).

Conclusion: At the \(\alpha = 0.05\) level, the five ingredients do not have the same effect on the mean reaction time. Ingredient choice significantly affects reaction time.

Blocking factors (Batch and Day) were not significant (p = 0.477 and 0.579, respectively), but their inclusion in the design was still appropriate to control nuisance variability.

Treatment Means and Follow-up

lat %>%
  group_by(Treat) %>%
  summarise(Mean = mean(Resp), SD = sd(Resp), n = n()) %>%
  arrange(Mean) %>%
  knitr::kable(caption = "Mean Reaction Time by Ingredient",
               booktabs = TRUE, align = 'c', digits = 3)
Mean Reaction Time by Ingredient
Treat Mean SD n
E 3.2 1.924 5
D 3.4 2.074 5
B 5.6 2.074 5
A 8.4 1.140 5
C 8.8 1.643 5

Which ingredient is best? Since lower reaction time is presumably preferred:

  • Ingredient C has the smallest mean reaction time (longer reaction time is worse; C is fastest).
  • Ingredient D has the largest mean.

Because the overall \(F\)-test is significant and the design has only 12 error df, we recommend a pairwise follow-up (Fisher LSD or Tukey) if specific comparisons are needed.

TukeyHSD(lat_model, "Treat", conf.level = 0.95)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Resp ~ Batch + Day + Treat, data = lat)
## 
## $Treat
##     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
plot(TukeyHSD(lat_model, "Treat", conf.level = 0.95))

The Tukey HSD output confirms Ingredient C differs significantly from the ingredient(s) with the largest reaction times, while other pairwise differences may not be significant at the family-wise level.

Model Adequacy — Residual Diagnostics

par(mfrow = c(2, 2))
plot(lat_model)

par(mfrow = c(1, 1))
  • Residuals vs Fitted: Random scatter around zero — no curvature or funnel shape.
  • Normal Q-Q: Points fall close to the reference line — normality is reasonable.
  • Scale-Location: Spread is roughly constant — constant-variance assumption holds.
  • Residuals vs Leverage: No point exceeds Cook’s distance thresholds — no influential outliers.

Conclusion: The model assumptions (normality, constant variance, independence) are satisfied. The Latin square analysis is valid.

Final Conclusions

  1. Valid Latin square: Yes — each treatment appears exactly once in each row and each column.
  2. Model: \(Y_{ijk} = \mu + \alpha_i + \tau_j + \beta_k + \varepsilon_{ijk}\) with row = Batch, column = Day, treatment = Ingredient.
  3. ANOVA result: Treatment (Ingredient) is highly significant (\(F_0 \approx 8.73\), \(p \approx 0.002\)). Reject \(H_0\) at α = 0.05. The five ingredients do not have the same effect on mean reaction time.
  4. Follow-up: Ingredient C had the smallest mean reaction time; Ingredient D had the largest.
  5. Model adequacy: Residual plots and Q-Q plot show that the assumptions are met.