# Analysis of variance table for 2x2 cross-over trials

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
close(con)

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
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
 396.8583


Also, the Treatment SS is calculated by

n1*n2/(2*(n1+n2))*(y11.bar-y12.bar-y21.bar+y22.bar)^2
 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?