Using Tukey's test for main effects to correct for multiplicity in a repeated measures ANOVA using R
(There is also example in the post-hoc talk demo file available at the Graduate statistics talk page from here)
install.packages(c("MASS")) install.packages(c("multcomp")) library(multcomp) library(MASS) library(car) library(nlme)
> oats
|
B |
V |
N |
Y |
|||
1 |
I |
Victory |
0.0cwt |
111 |
|||
2 |
I |
Victory |
0.2cwt |
130 |
|||
3 |
I |
Victory |
0.4cwt |
157 |
|||
4 |
I |
Victory |
0.6cwt |
174 |
|||
5 |
I |
Golden.rain |
0.0cwt |
117 |
|||
6 |
I |
Golden.rain |
0.2cwt |
114 |
|||
7 |
I |
Golden.rain |
0.4cwt |
161 |
|||
8 |
I |
Golden.rain |
0.6cwt |
141 |
|||
9 |
I |
Marvellous |
0.0cwt |
105 |
|||
10 |
I |
Marvellous |
0.2cwt |
140 |
|||
11 |
I |
Marvellous |
0.4cwt |
118 |
|||
12 |
I |
Marvellous |
0.6cwt |
156 |
|||
13 |
II |
Victory |
0.0cwt |
61 |
|||
14 |
II |
Victory |
0.2cwt |
91 |
|||
15 |
II |
Victory |
0.4cwt |
97 |
|||
16 |
II |
Victory |
0.6cwt |
100 |
|||
17 |
II |
Golden.rain |
0.0cwt |
70 |
|||
18 |
II |
Golden.rain |
0.2cwt |
108 |
|||
19 |
II |
Golden.rain |
0.4cwt |
126 |
|||
20 |
II |
Golden.rain |
0.6cwt |
149 |
|||
21 |
II |
Marvellous |
0.0cwt |
96 |
|||
22 |
II |
Marvellous |
0.2cwt |
124 |
|||
23 |
II |
Marvellous |
0.4cwt |
121 |
|||
24 |
II |
Marvellous |
0.6cwt |
144 |
|||
25 |
III |
Victory |
0.0cwt |
68 |
|||
26 |
III |
Victory |
0.2cwt |
64 |
|||
27 |
III |
Victory |
0.4cwt |
112 |
|||
28 |
III |
Victory |
0.6cwt |
86 |
|||
29 |
III |
Golden.rain |
0.0cwt |
60 |
|||
30 |
III |
Golden.rain |
0.2cwt |
102 |
|||
31 |
III |
Golden.rain |
0.4cwt |
89 |
|||
32 |
III |
Golden.rain |
0.6cwt |
96 |
|||
33 |
III |
Marvellous |
0.0cwt |
89 |
|||
34 |
III |
Marvellous |
0.2cwt |
129 |
|||
35 |
III |
Marvellous |
0.4cwt |
132 |
|||
36 |
III |
Marvellous |
0.6cwt |
124 |
|||
37 |
IV |
Victory |
0.0cwt |
74 |
|||
38 |
IV |
Victory |
0.2cwt |
89 |
|||
39 |
IV |
Victory |
0.4cwt |
81 |
|||
40 |
IV |
Victory |
0.6cwt |
122 |
|||
41 |
IV |
Golden.rain |
0.0cwt |
64 |
|||
42 |
IV |
Golden.rain |
0.2cwt |
103 |
|||
43 |
IV |
Golden.rain |
0.4cwt |
132 |
|||
44 |
IV |
Golden.rain |
0.6cwt |
133 |
|||
45 |
IV |
Marvellous |
0.0cwt |
70 |
|||
46 |
IV |
Marvellous |
0.2cwt |
89 |
|||
47 |
IV |
Marvellous |
0.4cwt |
104 |
|||
48 |
IV |
Marvellous |
0.6cwt |
117 |
|||
49 |
V |
Victory |
0.0cwt |
62 |
|||
50 |
V |
Victory |
0.2cwt |
90 |
|||
51 |
V |
Victory |
0.4cwt |
100 |
|||
52 |
V |
Victory |
0.6cwt |
116 |
|||
53 |
V |
Golden.rain |
0.0cwt |
80 |
|||
54 |
V |
Golden.rain |
0.2cwt |
82 |
|||
55 |
V |
Golden.rain |
0.4cwt |
94 |
|||
56 |
V |
Golden.rain |
0.6cwt |
126 |
|||
57 |
V |
Marvellous |
0.0cwt |
63 |
|||
58 |
V |
Marvellous |
0.2cwt |
70 |
|||
59 |
V |
Marvellous |
0.4cwt |
109 |
|||
60 |
V |
Marvellous |
0.6cwt |
99 |
|||
61 |
VI |
Victory |
0.0cwt |
53 |
|||
62 |
VI |
Victory |
0.2cwt |
74 |
|||
63 |
VI |
Victory |
0.4cwt |
118 |
|||
64 |
VI |
Victory |
0.6cwt |
113 |
|||
65 |
VI |
Golden.rain |
0.0cwt |
89 |
|||
66 |
VI |
Golden.rain |
0.2cwt |
82 |
|||
67 |
VI |
Golden.rain |
0.4cwt |
86 |
|||
68 |
VI |
Golden.rain |
0.6cwt |
104 |
|||
69 |
VI |
Marvellous |
0.0cwt |
97 |
|||
70 |
VI |
Marvellous |
0.2cwt |
99 |
|||
71 |
VI |
Marvellous |
0.4cwt |
119 |
|||
72 |
VI |
Marvellous |
0.6cwt |
121 |
Aov.mod <- aov(Y ~ N + V + Error(B/V), data = oats) Lme.mod <- lme(Y ~ N + V, random = ~1 | B/V, data = oats) summary(Aov.mod)
Error = B
|
Df |
Sum Sq |
Mean Sq |
F |
|||
Residuals |
5 |
15875 |
3175.1 |
|
Error = B:V
|
Df |
Sum Sq |
Mean Sq |
F |
|||
V |
2 |
1786.4 |
893.18 |
1.4853 |
|||
Residuals |
10 |
6013.3 |
601.33 |
|
Error = Within
|
Df |
Sum Sq |
Mean Sq |
F |
|||
N |
3 |
20020.5 |
6673.5 |
41.053 |
|||
Residuals |
51 |
8290.5 |
162.6 |
|
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
anova(Lme.mod)
|
numDF |
denDF |
F-value |
p-value |
||||
(Intercept) |
1 |
51 |
245.14566 |
<.0001 |
||||
N |
3 |
51 |
41.05276 |
<.0001 |
||||
V |
2 |
10 |
1.48534 |
0.2724 |
summary(Lme.mod)
Linear mixed-effects model fit by REML
- Data: oats
AIC |
BIC |
logLik |
||
586.0688 |
605.7756 |
-284.0344 |
Random effects:
- Formula: ~1 | B
- (Intercept)
StdDev: 14.64486
- Formula: ~1 | V %in% B
- (Intercept) Residual
StdDev: 10.47344 12.74987
Fixed effects: Y ~ N + V
|
Value |
Std.Error |
DF |
t-value |
p-value |
|||||
(Intercept) |
79.91667 |
8.220345 |
51 |
9.721814 |
0.0000 |
|||||
N0.2cwt |
19.50000 |
4.249956 |
51 |
4.588283 |
0.0000 |
|||||
N0.4cwt |
34.83333 |
4.249956 |
51 |
8.196164 |
0.0000 |
|||||
N0.6cwt |
44.00000 |
4.249956 |
51 |
10.353049 |
0.0000 |
|||||
VMarvellous |
5.29167 |
7.078905 |
10 |
0.747526 |
0.4720 |
|||||
VVictory |
-6.87500 |
7.078905 |
10 |
-0.971195 |
0.3544 |
Correlation:
(Intr) |
N0.2cw |
N0.4cw |
N06.cw |
VMrv11 |
|||
-0.259 |
- |
- |
- |
- |
|||
-0.259 |
0.500 |
- |
- |
- |
|||
-0.259 |
0.500 |
0.500 |
- |
- |
|||
-0.431 |
0.000 |
0.000 |
0.000 |
- |
|||
-0.431 |
0.000 |
0.000 |
0.000 |
0.500 |
Standardized Within-Group Residuals:
Min |
Q1 |
Med |
Q3 |
Max |
|||
-1.84134117 |
-0.66279719 |
-0.06694235 |
0.63822499 |
1.66 066729 |
Number of Observations: 72
Number of Groups:
|||||| B || V %in% B ||
6 |
18 |
summary(glht(Lme.mod, linfct=mcp(V="Tukey")))
- Simultaneous Tests for General Linear Hypotheses
Multiple Comparisons of Means: Tukey Contrasts
Fit: lme.formula(fixed = Y ~ N + V, data = oats, random = ~1 | B/V)
Linear Hypotheses:
|
Estimate |
Std. Error |
z value |
Pr(>|z|) |
||||
Marvellous - Golden.rain == 0 |
5.292 |
7.079 |
0.748 |
0.735 |
||||
Victory - Golden.rain == 0 |
-6.875 |
7.079 |
-0.971 |
0.595 |
||||
Victory - Marvellous == 0 |
-12.167 |
7.079 |
-1.719 |
0.198 |
(Adjusted p values reported -- single-step method)
Similarly for the within subjects effect
summary(glht(Lme.mod, linfct=mcp(N="Tukey")))