Latin Squares, part 1

An alternative to RCBDs that works with two nuisance factors.

Readings 4.2, 4.3, Rmarkdown

1. RCBD is useful when we have one nuisance factor. But what if we have two?

• We’re testing a manufacturing procedure, but raw materials come in batches and different operators have different skills.
• We’re testing different diets on cows over a series of days, but there will both cow and day effects.

2. Assume two nuisance factors, each with \(p\) levels. Furthermore, assume that we care about \(p\) different treatments.

3. The intuition for Latin Squares is similar to the intuition for RCBD, > Make sure to run each treatment on each nuisance factor exactly once.

   For example, in the manufacturing example, each operator should see each treatment once, and each batch of materials should be used for each treatment

Setup

4. Latin Squares are \(p\times p\) arrays, filled with \(p\) letters, so that each letter appears exactly once in each row and each column. For example, here is a table when \(p = 3\).

A	B	C
B	C	A
C	A	B

5. Why do we care? It tells us how we can implement the design idea above. First, randomly choose a \(p\times p\) Latin square
   • The rows are levels of the first blocking factor.
   • The columns are levels of the second blocking factor.
   • The letters are the treatment levels
   Then, the experiment consists of \(p^2\) runs, one for each of the pairs of blocking levels, with treatment specified by the cell’s label.

Model Setting

6. Instead of just one set of block effects, we’ll have two sets, \(\alpha_i\) and \(\beta_k\). This results in, \[\begin{align*} y_{ijk} &= \mu + \alpha_i + \tau_j + \beta_k + \epsilon_{ijk} \end{align*}\] where \(\epsilon_{ijk}\sim \mathcal{N}\left(0,\sigma^2\right)\) independently. Note that each of the indices \(i, j\) and \(k\) range from \(1, \dots, p\).

7. To make the model identifiable, we assume \[ \sum_{i = 1}^{p} \alpha_i = \sum_{j = 1}^{p} \tau_j = \sum_{k = 1}^{p} \beta_{k} = 0 \]

Hypothesis Testing

8. Our test hasn’t changed at all, \[\begin{align*} H_0 &: \tau_1 = \dots = \tau_a = 0 \\ H_{1} &: \tau_{i} \neq 0 \text{ for at least one } i \end{align*}\]

9. But now we need to account for block-to-block variation across both nuisance factors. The formula isn’t pretty, but it’s exactly analogous to the decompositions we’ve seen before,

\[ \begin{aligned} \sum_{i=1}^{p} \sum_{j=1}^{p} \sum_{k=1}^{p}\left(y_{i j k}-\bar{y} . . .\right)^{2}=& p \sum_{i=1}^{p}\left(\bar{y}_{i . .}-\bar{y} . . .\right)^{2}+\\ & p \sum_{j=1}^{p}\left(\bar{y}_{\cdot j} \cdot \bar{y}_{\cdots} .\right)^{2}+\\ & p \sum_{k=1}^{p}\left(\bar{y}_{\cdot \cdot k}-\bar{y} . .\right)^{2}+\\ & \sum_{i=1}^{p} \sum_{j=1}^{p} \sum_{k=1}^{p}\left(y_{i j k}-\bar{y}_{i . .}-\bar{y}_{. j .}+\bar{y} . .\right)^{2} \end{aligned} \]

10. This is abbreviated as, \[\begin{align*} SS_{\text{Total}} = &SS_{\text{Rows}} + \\ &SS_{\text{Columns}} +\\ &SS_{\text{Treatment}} + \\ &SS_{E} \end{align*}\] and we define

• \(MS_{\text{Rows}}=\frac{1}{p - 1}SS_{\text{Rows}}\)
• \(MS_{\text{Treatment}}=\frac{1}{p - 1}SS_{\text{Treatment}}\)
• \(MS_{\text{Columns}}=\frac{1}{p - 1}SS_{\text{Columns}}\)
• \(MS_{E}=\frac{1}{(p - 1)(p - 2)}SS_{E}\)

11. It turns out that \[ \frac{MS_{\text{Treatment}}}{MS_{E}} \sim F\left(p - 1, \left(p - 1\right)\left(p - 2\right)\right) \] which forms the basis for the test: reject the null if the ratio lies above the \(1 - \alpha\) quantile of this \(F\)-distribution.

Code Example

12. We’ll analyze the results of a study that used a latin square in its design. Compare the table below with table 4.9 in the book.

library(readr)
library(dplyr)
rocket <- read_table2("https://uwmadison.box.com/shared/static/ac68766l3zcjog1ldrzki3lis74bbd71.txt") %>%
  mutate_at(vars(-BurningRate), as.factor) # convert all but BurningRate to factor
rocket

# A tibble: 25 × 4
   Batch Operator Formulation BurningRate
   <fct> <fct>    <fct>             <dbl>
 1 1     1        a                    24
 2 1     2        b                    20
 3 1     3        c                    19
 4 1     4        d                    24
 5 1     5        e                    24
 6 2     1        b                    17
 7 2     2        c                    24
 8 2     3        d                    30
 9 2     4        e                    27
10 2     5        a                    36
# … with 15 more rows

13. Given this design, we can fit the model using a linear model. Here, \(\alpha_{i}, \tau_{j}\), and \(\beta_{k}\) are the batch, formulation, and operator, respectively. We’ll use the shorthand y ~ . to refer to the model using all the other variables in the data frame as inputs. Compare the ANOVA table below with Table 4.12.

#fit <- lm(BurningRate ~ Batch + Operator + Formulation, data = rocket) # gives exact same result
fit <- lm(BurningRate ~ ., data = rocket)
summary(aov(fit))

            Df Sum Sq Mean Sq F value  Pr(>F)   
Batch        4     68   17.00   1.594 0.23906   
Operator     4    150   37.50   3.516 0.04037 * 
Formulation  4    330   82.50   7.734 0.00254 **
Residuals   12    128   10.67                   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

13. There is an operator effect, but no detectable variation across batches. Controlling for batch and operater, there is a significant difference across formulations.