Latin Square Designs

Dynamite Example

> ### Latin Square Design, Dynamite Data
> ## Assume data are arranged as seen in the class notes
> x <- matrix(scan(what=""), 5, 10, byrow=TRUE)
1:  A 24   B 20   C 19   D 24   E 24
11:  B 17   C 24   D 30   E 27   A 36
21:  C 18   D 38   E 26   A 27   B 21
31:  D 26   E 31   A 26   B 23   C 22
41:  E 22   A 30   B 20   C 29   D 31
51: 
Read 50 items
> ( Latin <- x[,seq(1,9,2)] )
     [,1] [,2] [,3] [,4] [,5]
[1,] "A"  "B"  "C"  "D"  "E" 
[2,] "B"  "C"  "D"  "E"  "A" 
[3,] "C"  "D"  "E"  "A"  "B" 
[4,] "D"  "E"  "A"  "B"  "C" 
[5,] "E"  "A"  "B"  "C"  "D" 
> ( y <- matrix(as.numeric(x[,seq(2,10,2)]),5,5) )
     [,1] [,2] [,3] [,4] [,5]
[1,]   24   20   19   24   24
[2,]   17   24   30   27   36
[3,]   18   38   26   27   21
[4,]   26   31   26   23   22
[5,]   22   30   20   29   31
> ( batch <- row(y) )
     [,1] [,2] [,3] [,4] [,5]
[1,]    1    1    1    1    1
[2,]    2    2    2    2    2
[3,]    3    3    3    3    3
[4,]    4    4    4    4    4
[5,]    5    5    5    5    5
> ( operator <- col(y) )
     [,1] [,2] [,3] [,4] [,5]
[1,]    1    2    3    4    5
[2,]    1    2    3    4    5
[3,]    1    2    3    4    5
[4,]    1    2    3    4    5
[5,]    1    2    3    4    5
> dynamite <- data.frame(batch=factor(batch[1:25]),
+     operator=factor(operator[1:25]),
+     formula=factor(Latin[1:25]), response=y[1:25])
> head(dynamite)
  batch operator formula response
1     1        1       A       24
2     2        1       B       17
3     3        1       C       18
4     4        1       D       26
5     5        1       E       22
6     1        2       B       20
> tail(dynamite)
   batch operator formula response
20     5        4       C       29
21     1        5       E       24
22     2        5       A       36
23     3        5       B       21
24     4        5       C       22
25     5        5       D       31
> sapply(dynamite, class)
    batch  operator   formula  response 
 "factor"  "factor"  "factor" "numeric" 
> rm(x, y, Latin, batch, operator)   # tidy up
> save.image()   # save the image
> summary( dy.lm <- lm(response ~ ., data=dynamite) )

Call:
lm(formula = response ~ ., data = dynamite)

Residuals:
   Min     1Q Median     3Q    Max 
  -3.2   -1.2   -0.2    1.0    5.0 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   21.400      2.355   9.087 9.98e-07 ***
batch2         4.600      2.066   2.227  0.04586 *  
batch3         3.800      2.066   1.840  0.09067 .  
batch4         3.400      2.066   1.646  0.12568    
batch5         4.200      2.066   2.033  0.06475 .  
operator2      7.200      2.066   3.486  0.00450 ** 
operator3      2.800      2.066   1.356  0.20021    
operator4      4.600      2.066   2.227  0.04586 *  
operator5      5.400      2.066   2.614  0.02262 *  
formulaB      -8.400      2.066  -4.067  0.00156 ** 
formulaC      -6.200      2.066  -3.002  0.01103 *  
formulaD       1.200      2.066   0.581  0.57203    
formulaE      -2.600      2.066  -1.259  0.23207    
---
Signif. codes:  0 *** 0.001 ** 0.01 * 0.05 . 0.1   1

Residual standard error: 3.266 on 12 degrees of freedom
Multiple R-squared:  0.8107,    Adjusted R-squared:  0.6213 
F-statistic: 4.281 on 12 and 12 DF,  p-value: 0.008855

> anova(dy.lm)
Analysis of Variance Table

Response: response
          Df Sum Sq Mean Sq F value   Pr(>F)   
batch      4     68  17.000  1.5938 0.239059   
operator   4    150  37.500  3.5156 0.040373 * 
formula    4    330  82.500  7.7344 0.002537 **
Residuals 12    128  10.667                    
---
Signif. codes:  0 *** 0.001 ** 0.01 * 0.05 . 0.1   1
> ( dy.tk <- TukeyHSD(aov(response ~ ., data=dynamite), "formula") )
  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = response ~ ., data = dynamite)

$formula
    diff         lwr        upr     p adj
B-A -8.4 -14.9839317 -1.8160683 0.0110827
C-A -6.2 -12.7839317  0.3839317 0.0684350
D-A  1.2  -5.3839317  7.7839317 0.9754380
E-A -2.6  -9.1839317  3.9839317 0.7194121
C-B  2.2  -4.3839317  8.7839317 0.8204614
D-B  9.6   3.0160683 16.1839317 0.0041583
E-B  5.8  -0.7839317 12.3839317 0.0944061
D-C  7.4   0.8160683 13.9839317 0.0254304
E-C  3.6  -2.9839317 10.1839317 0.4461852
E-D -3.8 -10.3839317  2.7839317 0.3966727

> windows(width=5, height=5, pointsize=10)
> plot(dy.tk, las=1)
> mtext("Tukey's Honest Significant Difference", side=3, line=0,
+       col="NavyBlue")

The Tukey's HSD multiple comparison results are graphed:
(dynamite01.png here)