Chapter 41 Latin Squares
41.1 Data and Preliminaries
Data from Kirk (1982, Brooks/Cole). Adjective/noun pairs are presented on a computer screen to subjects. The adjective is present for 100ms, followed by the noun for 20 ms. If the subject does not recognize the noun, the pair is presented again with the presentation time for the noun increased by 5 ms. This is repeated until the noun is recognized. Of interest is whether the degree of association between the adjective and the noun influences how quickly the subject recognizes the noun. For example, you might recognize “ice cream” more quickly after “vanilla” than after “purple”.
Twenty-five subjects are used, and the subjects differ in their overall quickness to recognize nouns, so the subjects are blocked into five blocks of size five based on how quick they are. The experiment will be administered by five graduate students, who may have some effect, so the subjects are also blocked by which student administers the test. The response measured is the average time in ms taken to recognize 10 nouns of a given association strength. They used a Latin square because of the two blocking factors.
> library(cfcdae)
Registered S3 method overwritten by 'DoE.base':
method from
factorize.factor conf.design
> library(lme4)
Loading required package: Matrix
> library(car)
Loading required package: carData
> nouns <- read.table("nouns.dat.txt",header=TRUE)
> nouns <- within(nouns,{oquick<-factor(oquick);student<-factor(student);assocstr<-factor(assocstr)})
> nouns
oquick student assocstr rtime
1 1 1 1 72
2 1 2 2 62
3 1 3 3 66
4 1 4 4 51
5 1 5 5 40
6 2 1 2 65
7 2 2 3 61
8 2 3 4 40
9 2 4 5 44
10 2 5 1 59
11 3 1 3 55
12 3 2 4 46
13 3 3 5 35
14 3 4 1 63
15 3 5 2 54
16 4 1 4 34
17 4 2 5 29
18 4 3 1 54
19 4 4 2 44
20 4 5 3 50
21 5 1 5 51
22 5 2 1 49
23 5 3 2 43
24 5 4 3 30
25 5 5 4 25
Here is the acutal LS used; it is perhaps the most boring LS possible. Rows are blocks for overall quickness, columns are blocks for students, and matrix elements are the treatments (strength of association between adjective and noun). There’s likely a cleaner, prettier way to get this table, but this works.
> with(nouns,tapply(as.numeric(assocstr),list(oquick,student),mean))
1 2 3 4 5
1 1 2 3 4 5
2 2 3 4 5 1
3 3 4 5 1 2
4 4 5 1 2 3
5 5 1 2 3 4
41.1.1 Fixed Blocks
Fit using fixed blocks and plot residuals. It’s not the best looking plot. Variability decreasing with mean? Maybe a little curvature?
> nouns.lm <- lm(rtime~oquick+student+assocstr,data=nouns)
> plot(nouns.lm,which=1)

Figure 41.1: Residual plot of nouns LS
Normal plot looks slightly long tailed. Just a hint of evidence that #21 is an outlier.
> plot(nouns.lm,which=2)

Figure 41.2: NPP of residuals for nouns LS
Let’s follow up on the issue of outliers using the outlier test. There is some evidence, but not overwhelming evidence, that #21 does not follow the pattern of the rest of the data.
> outlierTest(nouns.lm)
rstudent unadjusted p-value Bonferroni p
21 4.23749 0.0013947 0.034869
One thing that we can do with an outlier is to make a dummy variable that points just to the outlier, and then include the dummy variable in the model. Then the other data get fit as if the outlier weren’t even in the data set, and the coefficient of the dummy variable tells us how far the outlier is from what would be predicted by the remainder of the data. This is a slightly more informative approach than simply deleting an outlier.
Update the model to
include the dummy and look at what we got. This indicates that the
outlier is about 29 ms higher than the rest of the data would have
predicted. Note that the t-value for odd
matches the studentized residual
we had for response 21 in the previous model.
> odd <- rep(0,25);odd[21] <- 1
> nouns.lm.odd <- update(nouns.lm,~odd+.)
> summary(nouns.lm.odd)
Call:
lm(formula = rtime ~ odd + oquick + student + assocstr, data = nouns)
Residuals:
Min 1Q Median 3Q Max
-7.1500 -2.1500 0.2333 2.2333 5.2500
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 47.7167 0.9898 48.206 3.75e-14 ***
odd 29.0833 6.8633 4.237 0.001395 **
oquick1 10.4833 1.9217 5.455 0.000199 ***
oquick2 6.0833 1.9217 3.166 0.008990 **
oquick3 2.8833 1.9217 1.500 0.161661
oquick4 -5.5167 1.9217 -2.871 0.015223 *
student1 1.8667 2.1963 0.850 0.413483
student2 1.6833 1.9217 0.876 0.399784
student3 -0.1167 1.9217 -0.061 0.952680
student4 -1.3167 1.9217 -0.685 0.507433
assocstr1 11.6833 1.9217 6.080 7.97e-05 ***
assocstr2 5.8833 1.9217 3.061 0.010824 *
assocstr3 4.6833 1.9217 2.437 0.032996 *
assocstr4 -8.5167 1.9217 -4.432 0.001009 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 4.755 on 11 degrees of freedom
Multiple R-squared: 0.9344, Adjusted R-squared: 0.8569
F-statistic: 12.06 on 13 and 11 DF, p-value: 0.0001094
Recheck residuals. These look better than before: less “humped in the middle” appearance and seemingly more constant variance.
> plot(nouns.lm.odd,which=1)

Figure 41.3: Residual plot of nouns data with #21 fit separately
Normality looks better too.
> plot(nouns.lm.odd,which=2)
Warning: not plotting observations with leverage one:
21

Figure 41.4: NPP of residuals in nouns data with #21 fit separately
At this point all we can say is the #21 does not fit the pattern of the rest of the data. We cannot say that it is a bad data point. We need to go back and consult the lab notebooks to see if something went wrong.
We will continue fitting with #21 identified separately.
Things are unbalanced and messier when we remove outliers or fit them individually. We can use Type II SS. The result here for treatment is the same as in the preceding sequential anova, because treatment was entered last.
> Anova(nouns.lm.odd,type=2)
Anova Table (Type II tests)
Response: rtime
Sum Sq Df F value Pr(>F)
odd 406.00 1 17.9563 0.001395 **
oquick 1595.62 4 17.6424 9.431e-05 ***
student 57.28 4 0.6334 0.649113
assocstr 1976.23 4 21.8507 3.435e-05 ***
Residuals 248.72 11
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
The strongest associations (levels 4 and 5) cannot be distinguished from each other, and the weakest (levels 1, 2, and 3) cannot be distinguished from each other, but the strongest can be distinguished from the weakest.
> sidelines(pairwise(nouns.lm.odd,assocstr))
5 -13.73 |
4 -8.52 |
3 4.68 |
2 5.88 |
1 11.68 |
Let’s look at the treatment effects. They are monotone decreasing.
> model.effects(nouns.lm.odd,"assocstr")
1 2 3 4 5
11.683333 5.883333 4.683333 -8.516667 -13.733333
We don’t know that these strengths in any sense equally spaced. However, if we assume for the moment that they are, we could plot the treatment effects against 1 through 5; Except for level 3, this looks incredibly linear (too linear to be real?).
> plot(1:5,model.effects(nouns.lm.odd,"assocstr"))

Figure 41.5: Plot of association strength affects against level
Linear, quadratic, cubic, and quartic contrasts for equally spaced, equal sample size treatments. These are not quite correct here, because we have fit #21 separately.
Linear is significant; nothing else is. The joint
(combined) F is the same as we saw for
assocstr
in the Type II ANOVA above.
> cfs <- matrix(c(-2,-1,0,1,2, 2,-1,-2,-1,2, -1,2,0,-2,1, 1,-4,6,-4,1),nrow=5)
> cfs
[,1] [,2] [,3] [,4]
[1,] -2 2 -1 1
[2,] -1 -1 2 -4
[3,] 0 -2 0 6
[4,] 1 -1 -2 -4
[5,] 2 2 1 1
> linear.contrast(nouns.lm.odd,assocstr,cfs,joint=TRUE)
$estimates
estimates se t-value p-value lower-ci upper-ci
1 -65.233333 7.263476 -8.9810079 2.139925e-06 -81.220136 -49.246531
2 -10.833333 8.417040 -1.2870716 2.244931e-01 -29.359114 7.692447
3 3.383333 6.863340 0.4929573 6.317354e-01 -11.722775 18.489442
4 36.583333 17.844683 2.0500972 6.496480e-02 -2.692549 75.859215
$Ftest
F df1 df2 p-value
21.85073 4 11 3.434593e-05
attr(,"class")
[1] "linear.contrast"
If we want to test whether there is anything beyond linear in strength (and do it right), then fit a model with linear and compare it to the full model. There is no evidence of anything beyond linear being needed.
> nouns.lm.odd.linear <- lm(rtime~odd+oquick+student+as.numeric(assocstr),data=nouns)
> anova(nouns.lm.odd.linear,nouns.lm.odd)
Analysis of Variance Table
Model 1: rtime ~ odd + oquick + student + as.numeric(assocstr)
Model 2: rtime ~ odd + oquick + student + assocstr
Res.Df RSS Df Sum of Sq F Pr(>F)
1 14 391.11
2 11 248.72 3 142.39 2.0992 0.1585
41.1.2 Random blocks
One might hypothesize that the graduate students used
to administer the test are just a random
sample of potential graduate students, thereby indicating
that student
should be random.
If you fit with random student
without separating #21, you get
the same thing as we got for fixed student
. (Anova shown,
summaries would be alike as well.)
> nouns.lmer <- lmer(rtime~oquick+(1|student)+assocstr,data=nouns)
> Anova(nouns.lmer,test="F")
Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
Response: rtime
F Df Df.res Pr(>F)
oquick 5.6078 4 12 0.008798 **
assocstr 7.3655 4 12 0.003094 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> anova(nouns.lm)
Analysis of Variance Table
Response: rtime
Df Sum Sq Mean Sq F value Pr(>F)
oquick 4 1223.84 305.96 5.6078 0.008798 **
student 4 306.64 76.66 1.4051 0.290720
assocstr 4 1607.44 401.86 7.3655 0.003094 **
Residuals 12 654.72 54.56
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
However, if you fit #21 separately, then the random
and fixed models give similar, but not identical, results.
The loss of balance/orthogonality leads to differences.
Here we see that the F for assocstr
(in the model with odd
)
is larger under lmer()
than under lm()
.
> nouns.lmer.odd <- lmer(rtime~odd+oquick+(1|student)+assocstr,data=nouns)
boundary (singular) fit: see help('isSingular')
> Anova(nouns.lmer.odd,test="F",type=2)
Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
Response: rtime
F Df Df.res Pr(>F)
odd 24.986 1 15.00 0.0001588 ***
oquick 20.418 4 11.33 3.954e-05 ***
assocstr 25.067 4 11.33 1.435e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> Anova(nouns.lm.odd,type=2)
Anova Table (Type II tests)
Response: rtime
Sum Sq Df F value Pr(>F)
odd 406.00 1 17.9563 0.001395 **
oquick 1595.62 4 17.6424 9.431e-05 ***
student 57.28 4 0.6334 0.649113
assocstr 1976.23 4 21.8507 3.435e-05 ***
Residuals 248.72 11
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
You can see an interesting thing when you compare the estimated
variances in the lmer()
models with and without odd
.
Without odd
, we estimate a positive variance for student
and \(\hat{\sigma}^2 = 54.56\). If you fit odd
, then lmer()
estimates
a zero variance for student
and \(\hat{\sigma}^2 = 20.4\).
The estimated variance for lm()
using odd
was 22.61.
> VarCorr(nouns.lmer)
Groups Name Std.Dev.
student (Intercept) 2.1024
Residual 7.3865
> VarCorr(nouns.lmer.odd)
Groups Name Std.Dev.
student (Intercept) 0.0000
Residual 4.5166
lmer()
thinks that student
has 0 variability when odd
is in the model. But there is actually some variability
between students, just not enough to give a positive REML
estimate. If you combine the SS and DF for student
and
residual from the lm()
model with odd
, you get
something that closely matches the lmer()
residual
variance in the model with odd
.
> sqrt((248.72+57.28)/15)
[1] 4.516636
41.2 Milk data
These data give the milk production (in pounds) for 18 cows for three consecutive five-week periods. Each cow is given a different diet in each period. The diets are good hay, poor hay, and straw. Data from John (1971 Wiley).
The experiment is arranged as a cross-over design – six Latin squares with cows as column blocks and time periods as row blocks.
> milkdata <- read.table("milk.dat.txt",header=TRUE)
> milkdata <- within(milkdata,{cow <- factor(cow);allcows <- factor(allcows);
+ period <- factor(period);square <- factor(square);diet <- factor(diet)})
> milkdata
cow period square diet allcows milk
1 1 1 1 1 1 768
2 1 2 1 2 1 600
3 1 3 1 3 1 411
4 2 1 1 2 2 662
5 2 2 1 3 2 515
6 2 3 1 1 2 506
7 3 1 1 3 3 731
8 3 2 1 1 3 680
9 3 3 1 2 3 525
10 1 1 2 1 4 669
11 1 2 2 3 4 550
12 1 3 2 2 4 416
13 2 1 2 2 5 459
14 2 2 2 1 5 409
15 2 3 2 3 5 222
16 3 1 2 3 6 624
17 3 2 2 2 6 462
18 3 3 2 1 6 426
19 1 1 3 1 7 1091
20 1 2 3 2 7 798
21 1 3 3 3 7 534
22 2 1 3 2 8 1234
23 2 2 3 3 8 902
24 2 3 3 1 8 869
25 3 1 3 3 9 1300
26 3 2 3 1 9 1297
27 3 3 3 2 9 962
28 1 1 4 1 10 1105
29 1 2 4 3 10 712
30 1 3 4 2 10 453
31 2 1 4 2 11 891
32 2 2 4 1 11 830
33 2 3 4 3 11 629
34 3 1 4 3 12 859
35 3 2 4 2 12 617
36 3 3 4 1 12 597
37 1 1 5 1 13 941
38 1 2 5 2 13 718
39 1 3 5 3 13 548
40 2 1 5 2 14 794
41 2 2 5 3 14 603
42 2 3 5 1 14 613
43 3 1 5 3 15 779
44 3 2 5 1 15 718
45 3 3 5 2 15 515
46 1 1 6 1 16 933
47 1 2 6 3 16 658
48 1 3 6 2 16 576
49 2 1 6 2 17 724
50 2 2 6 1 17 649
51 2 3 6 3 17 496
52 3 1 6 3 18 749
53 3 2 6 2 18 594
54 3 3 6 1 18 612
41.2.1 Fixed effects
The simple approach is to just use lm()
to fit the Latin
Square. We have two choices for how to model the cows. We
can enumerate the cows 1 through 18 (allcows
) or we can
use square and cow within square (square/cow
). Note that
the sum of SS for square
and cow:square
is the same
as that for allcows
(df add as well).
> milk.lm.squares <- lm(milk~square/cow+period+diet,data=milkdata)
> milk.lm.allcows <- lm(milk~allcows+period+diet,data=milkdata)
> anova(milk.lm.squares)
Analysis of Variance Table
Response: milk
Df Sum Sq Mean Sq F value Pr(>F)
square 5 1392534 278507 80.6099 < 2.2e-16 ***
period 2 814222 407111 117.8326 1.742e-15 ***
diet 2 121147 60573 17.5321 7.220e-06 ***
square:cow 12 318241 26520 7.6759 1.941e-06 ***
Residuals 32 110560 3455
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> anova(milk.lm.allcows)
Analysis of Variance Table
Response: milk
Df Sum Sq Mean Sq F value Pr(>F)
allcows 17 1710775 100634 29.127 9.384e-15 ***
period 2 814222 407111 117.833 1.742e-15 ***
diet 2 121147 60573 17.532 7.220e-06 ***
Residuals 32 110560 3455
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> 1392534+318241
[1] 1710775
Checking residuals, we see the flopping fish: curvature and non-constant variance. We need a transformation.
> plot(milk.lm.squares,which=1)

Figure 41.6: Residual plot of milk data
Box-Cox suggests something around .3; .5 is barely in the interval and 0 (log) is well outside the interval.
> boxCox(milk.lm.squares)

Figure 41.7: Box-Cox plot for milk data
We’ll try power .25. Residuals are better, but not great, with perhaps an outlier, which is confirmed by the outlier test.
> milk.lm.squares.25 <- lm(milk^.25 ~ square/cow+period+diet,data=milkdata)
> plot(milk.lm.squares.25,which=1)

Figure 41.8: Residual plot of fourth root milk data
> outlierTest(milk.lm.squares.25)
rstudent unadjusted p-value Bonferroni p
30 -4.134151 0.00025126 0.013568
Point #30 was the potential outlier, so create a dummy variable, refit, and recheck residuals. The residual plot is a little better, although there is one seemingly isolated point, and the outlier test does not find additional outliers.
> odd30 <- rep(0,54);odd30[30] <- 1
> milk.lm.squares.25.odd <- update(milk.lm.squares.25,~odd30+.)
> plot(milk.lm.squares.25.odd,which=1)

Figure 41.9: Residual plot for fourth power milk data with outlier
> outlierTest(milk.lm.squares.25.odd)
No Studentized residuals with Bonferroni p < 0.05
Largest |rstudent|:
rstudent unadjusted p-value Bonferroni p
21 -2.422427 0.021667 NA
OK, this model seems to fit well enough. Diet 1 is better than diet 2, and diet 3 must be worse still. Milk production is consistently decreasing over time, so that was a useful block.
The outlier is about .41 low (on the quarter power scale), which is a long way from where it should be (bigger than treatment effects, period effects, and nearly all the cow effects). Number 30 is cow 10, period 3, treatment 2. We need to check this to see if there is some reason that it is different.
> summary(milk.lm.squares.25.odd)
Call:
lm(formula = milk^0.25 ~ odd30 + square + period + diet + square:cow,
data = milkdata)
Residuals:
Min 1Q Median 3Q Max
-0.132524 -0.035098 -0.005765 0.049114 0.124407
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.09285 0.01057 481.862 < 2e-16 ***
odd30 -0.41073 0.09935 -4.134 0.000251 ***
square1 -0.16084 0.02335 -6.890 1.00e-07 ***
square2 -0.46977 0.02335 -20.123 < 2e-16 ***
square3 0.49391 0.02335 21.157 < 2e-16 ***
square4 0.14436 0.02502 5.769 2.37e-06 ***
square5 0.02005 0.02335 0.859 0.397025
period1 0.27635 0.01483 18.630 < 2e-16 ***
period2 -0.01140 0.01483 -0.769 0.447814
diet1 0.11449 0.01483 7.718 1.05e-08 ***
diet2 -0.02932 0.01517 -1.933 0.062476 .
square1:cow1 -0.02665 0.03605 -0.739 0.465402
square2:cow1 0.19182 0.03605 5.320 8.57e-06 ***
square3:cow1 -0.29700 0.03605 -8.238 2.64e-09 ***
square4:cow1 0.08122 0.04228 1.921 0.063947 .
square5:cow1 0.07154 0.03605 1.984 0.056135 .
square6:cow1 0.09835 0.03605 2.728 0.010403 *
square1:cow2 -0.07234 0.03605 -2.006 0.053601 .
square2:cow2 -0.29450 0.03605 -8.169 3.17e-09 ***
square3:cow2 0.02544 0.03605 0.706 0.485699
square4:cow2 0.04243 0.03771 1.125 0.269133
square5:cow2 -0.03306 0.03605 -0.917 0.366233
square6:cow2 -0.08055 0.03605 -2.234 0.032815 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.07648 on 31 degrees of freedom
Multiple R-squared: 0.9801, Adjusted R-squared: 0.966
F-statistic: 69.46 on 22 and 31 DF, p-value: < 2.2e-16
Treatment 1 (good hay) is different from the other two, but they cannot be distinguished (just barely at the .05 level with HSD).
> pairwise(milk.lm.squares.25.odd,diet)
Pairwise comparisons ( hsd ) of diet
estimate signif diff lower upper
* 1 - 2 0.14380950 0.06419838 0.079611120 0.2080079
* 1 - 3 0.19965943 0.06274465 0.136914789 0.2624041
2 - 3 0.05584993 0.06419838 -0.008348451 0.1200483
41.2.2 Random effects
What if we thought the cows were a random sample? Then
we would have the cow effect as a random effect.
For fixed effects, whether you used allcows
or
square/cow
made no difference; it makes a very
slight difference for random effects.
> milk.lmer.squares.25.odd <- lmer(milk^.25~period+diet+odd30+(1|square/cow),data=milkdata)
> milk.lmer.allcows.25.odd <- lmer(milk^.25~period+diet+odd30+(1|square:cow),data=milkdata)
> Anova(milk.lmer.squares.25.odd,test="F",type=2)
Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
Response: milk^0.25
F Df Df.res Pr(>F)
period 216.835 2 31.045 < 2.2e-16 ***
diet 32.600 2 31.045 2.358e-08 ***
odd30 16.482 1 32.543 0.0002893 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> Anova(milk.lmer.allcows.25.odd,test="F",type=2)
Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
Response: milk^0.25
F Df Df.res Pr(>F)
period 216.549 2 31.017 < 2.2e-16 ***
diet 32.572 2 31.017 2.394e-08 ***
odd30 16.648 1 31.573 0.0002847 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
The variance components are quite different in the two approaches. Using squares, there is large square to square variability and smaller cow within square variability. If these were just 18 cows randomly arranged into squares, we would not expect that. If these were six squares representing three cows at each of six farms, then it might not be surprising. Using just cows, the cow to cow variability has to be a lot bigger to account for variability between squares.
> VarCorr(milk.lmer.squares.25.odd)
Groups Name Std.Dev.
cow:square (Intercept) 0.167595
square (Intercept) 0.304048
Residual 0.076463
> VarCorr(milk.lmer.allcows.25.odd)
Groups Name Std.Dev.
square:cow (Intercept) 0.331195
Residual 0.076475
In the context of random, hierarchical blocks, it makes some sense to test whether higher levels of the hierarchy contribute variance. The square effect is rather significant.
> milk.lmer.onlysquares.25.odd <- lmer(milk^.25~period+diet+odd30+(1|square),data=milkdata)
> library(RLRsim)
> exactRLRT(milk.lmer.onlysquares.25.odd,milk.lmer.squares.25.odd,milk.lmer.allcows.25.odd)
simulated finite sample distribution of RLRT.
(p-value based on 10000 simulated values)
data:
RLRT = 10.684, p-value = 5e-04
Here we compare the models of “only square effects”, “only cow effects”, and “both square and cow effects” via AIC. The model with both levels of variation is clearly chosen, in agreement with the test we just ran.
> AIC(milk.lmer.onlysquares.25.odd)
[1] 17.33312
> AIC(milk.lmer.allcows.25.odd)
[1] -8.533075
> AIC(milk.lmer.squares.25.odd)
[1] -17.2169
41.3 Relative Efficiency
Let’s go back to the association strength data. Repeat the ANOVA.
> Anova(nouns.lm.odd,type=2)
Anova Table (Type II tests)
Response: rtime
Sum Sq Df F value Pr(>F)
odd 406.00 1 17.9563 0.001395 **
oquick 1595.62 4 17.6424 9.431e-05 ***
student 57.28 4 0.6334 0.649113
assocstr 1976.23 4 21.8507 3.435e-05 ***
Residuals 248.72 11
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
What if we had not blocked on original quickness? How much larger of an experiment would we have needed?
The MSE in the Latin Square is \(248.72/11 = 22.61\). The estimated MSE for an RCB that used only student is a weighted average of MSE and MS for original quickness. Note here that I have used the df for a square without an outlier being fitted. (I also note the futility of multiplying by four and then dividing by 4, but the formulae in the book and notes are based on mean squares.)
Clearly, original quickness made a big difference.
> 248.72/11
[1] 22.61091
> ((12+4)*(248.72/11)+4*(1595.62/4))/(12+4+4)
[1] 97.86973
Here is the relative efficiency. We would have needed four times as much data if we had not blocked on original quickness.
> (13/15)*(19/17)*97.871/22.61
[1] 4.192859
If you repeat this for student
, you will find
no gain.
Now look at the milk data. What is the relative efficiency if we had blocked on cow but not on period?
First estimate the variance from the ANOVA output.
> Anova(milk.lm.squares.25.odd,type=2)
Anova Table (Type II tests)
Response: milk^0.25
Sum Sq Df F value Pr(>F)
odd30 0.1000 1 17.091 0.0002513 ***
square 4.5647 5 156.076 < 2.2e-16 ***
period 2.5290 2 216.181 < 2.2e-16 ***
diet 0.3806 2 32.535 2.432e-08 ***
square:cow 1.0727 12 15.283 8.242e-10 ***
Residuals 0.1813 31
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> ((31+2)*(.1813/31)+2.5290)/(31+2+2)
[1] 0.07777134
This error is 13 times the error variance when we included period.
Now the final relative efficiency. There are 32 df for this replicated Latin square and 34 for the corresponding RCB, so the df adjustment does very little. The relative efficiency is over 13.
> (32+1)/(32+3) * (34+3)/(34+1)
[1] 0.9967347
> (32+1)/(32+3) * (34+3)/(34+1) * .0778 / (.1813/31)
[1] 13.25938