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)
Residual plot of nouns LS

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)
NPP of residuals for nouns LS

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)
Residual plot of nouns data with #21 fit separately

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
NPP of residuals in nouns data with #21 fit separately

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"))
Plot of association strength affects against level

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)
Residual plot of milk data

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)
Box-Cox plot for milk data

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)
Residual plot of fourth root milk data

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)
Residual plot for fourth power milk data with outlier

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