Analysis of variance table for 2x2 cross-over trials

Question

I am trying to reproduce the anova table for a 2x2 cross-over design. I used the data listed in tables 2.1 and 2.2 of the "Design and Analysis of cross-over trials" book, by Jones and Kenward. I am giving all the code so that you can reproduce the calculations.

g1AB <- read.table(textConnection("
Label Per1 Per2
7 121.905 116.667
8 218.5 200.5
9 235 217.143
13 250 196.429
14 186.19 185.5
15 231.563 221.842
17 443.25 420.5
21 198.421 207.692
22 270.5 213.158
28 360.476 384
35 229.75 188.25
36 159.091 221.905
37 255.882 253.571
38 279.048 267.619
41 160.556 163
44 172.105 182.381
58 267 313
66 230.75 211.111
71 271.19 257.619
76 276.25 222.105
79 398.75 404
80 67.778 70.278
81 195 223.158
82 325 306.667
86 368.077 362.5
89 228.947 227.895
90 236.667 220
")->con,header=T)
close(con)

g2BA <- read.table(textConnection("
Label Per1 Per2
3 138.333 138.571
10 225 256.25
11 392.857 381.429
16 190 233.333
18 191.429 228
23 226.19 267.143
24 201.905 193.5
26 134.286 128.947
27 238 248.5
29 159.5 140
30 232.75 276.563
32 172.308 170
33 266 305
39 171.333 186.333
43 194.737 191.429
47 200 222.619
51 146.667 183.81
52 208 241.667
55 208.75 218.81
59 271.429 225
68 143.81 188.5
70 104.444 135.238
74 145.238 152.857
77 215.385 240.476
78 306 288.333
83 160.526 150.476
84 353.81 369.048
85 293.889 308.095
99 371.190 404.762
")->con,header=T)
close(con)

n1 <- nrow(g1AB)
n2 <- nrow(g2BA)
p <- 2 # periods

Some quantities that may be useful:

y11. <- sum(g1AB$Per1)
y12. <- sum(g1AB$Per2)
y21. <- sum(g2BA$Per1)
y22. <- sum(g2BA$Per2)
y1.. <- sum(c(g1AB$Per1,g1AB$Per2))
y2.. <- sum(c(g2BA$Per1,g2BA$Per2))
y... <- y1.. + y2..

y11.bar <- 1/n1 * y11.
y12.bar <- 1/n1 * y12.
y21.bar <- 1/n2 * y21.
y22.bar <- 1/n2 * y22.

y1..bar <- 1/(p*n1)*(y11.+y12.)
y2..bar <- 1/(p*n2)*(y21.+y22.)

y...bar <- 1/(p*(n1+n2)) * (y1.. + y2..)

In order to perform the analysis of variance in R, I created the following dataset

mydata1 <- data.frame(PEFR=g1AB$Per1,Subjects=g1AB$Label,Time=1,Groups="AB",Treatment=1)
mydata2 <- data.frame(PEFR=g2BA$Per1,Subjects=g2BA$Label,Time=1,Groups="BA",Treatment=2)
mydata3 <- data.frame(PEFR=g1AB$Per2,Subjects=g1AB$Label,Time=2,Groups="AB",Treatment=2)
mydata4 <- data.frame(PEFR=g2BA$Per2,Subjects=g2BA$Label,Time=2,Groups="BA",Treatment=1)
mydata <- rbind(mydata1,mydata2,mydata3,mydata4)
mydata$Subjects<-factor(mydata$Subjects)

(I have read that) The correct anova table is obtained using the following command

res <- summary(PEFR.aov <- aov(PEFR~Groups+Time+Treatment+Error(Subjects), data=mydata))

Error: Subjects
          Df Sum Sq Mean Sq F value Pr(>F)
Groups     1  10573   10573   0.899  0.347
Residuals 54 634866   11757               

Error: Within
          Df Sum Sq Mean Sq F value  Pr(>F)   
Time       1    480   479.6   1.470 0.23061   
Treatment  1   3026  3026.1   9.276 0.00359 **
Residuals 54  17617   326.2 

1) Can you please explain to me the use of the Error term in the formula?

2) The anova table is the same to that of table 2.10 of the book, except the SS for Time (Period) which has a value of 396.858 instead of 480. According to Table 2.8 this quantity is calculated using

> n1*n2/(2*(n1+n2))*(y11.bar-y12.bar+y21.bar-y22.bar)^2
[1] 396.8583

Also, the Treatment SS is calculated by

n1*n2/(2*(n1+n2))*(y11.bar-y12.bar-y21.bar+y22.bar)^2
[1] 3026.12

If the formula was specified like this

PEFR~Groups+Treatment+Time+Error(Subjects)

the SS for Time would be correct, but the Treatment SS would have a value of 3109. How can I get 3026 for Treatment and 397 for Time, as I see in the book? What should be taken care of in the order the variables are specified in the formula?

Thank you for your time.

Answer 1

1) The Error term basically tells the aov() function that in addition to randomness between each observation, some of the observations are related to eachother by being observations on the same subject. So a random component of the PEFR value comes from each subject, before calculating the random component from each individual observation on that subject.

2) Your issue with the sums of squares looks like it is related to the question of Type III sums of squares (the default in some other stats packages). Type III SS are the marginal sum of squares explained by one variable or group of variables after all the variance explained by the other variables taken into account. As you have discovered, R only reports this for the last variable entered into the model. Although sometimes controversial, there are good reasons for this and I myself prefer the R default - I like my sums of squares in an ANOVA table to add up to the right number. As you have discovered, there are plenty of ways to generate the marginal sum of squares of a group of variables if you want to test them.

To answer your question what care should be taken in the order... you should put the variables you are interested in testing last. That way you are implicitly testing two nested models; one with all the variables including the ones in question; another with only the various other variables (eg nuisance confounding variables, or ones you are just positive are in the best model and hence don't need to be tested). You can then conduct a hypothesis test on the sum of squares reported for the last variable/s.