I'm working on mixed design ANOVA and would like to run TukeyHSD for its post-Hoc test.
I keep getting error, "Error in UseMethod("TukeyHSD") :
no applicable method for 'TukeyHSD' applied to an object of class "c('aovlist', 'listof')".
I've seen several questions about this errors in some F&Q websites, but I still cannot find a solution.
My dataset looks like this:
subject response time group
1 S1 0.99676 time1 task1
2 S2 1.00220 time1 task1
3 S3 1.00420 time1 task1
4 S4 0.99467 time1 task1
5 S5 0.97906 time1 task1
6 S6 0.99162 time1 task1
7 S7 0.99771 time1 task1
8 S8 1.01780 time1 task2
9 S9 0.98682 time1 task2
10 S10 0.99124 time1 task2
11 S11 1.01670 time1 task2
12 S12 0.99769 time1 task2
13 S13 1.02090 time1 task2
14 S14 1.01740 time1 task2
15 S15 0.98851 time1 task3
16 S16 1.00690 time1 task3
17 S17 0.99717 time1 task3
18 S18 0.98945 time1 task3
19 S19 1.00270 time1 task3
20 S20 1.02690 time1 task3
21 S21 1.00050 time1 task3
22 S1 0.96908 time2 task1
23 S2 0.94024 time2 task1
......
Then, the anova result is:
anova = aov(response~(group*time)+Error(subject/time),data) summary(anova)
Error: subject
Df Sum Sq Mean Sq F value Pr(>F)
group 2 0.00381 0.001907 0.701 0.509
Residuals 18 0.04896 0.002720
Error: subject:time
Df Sum Sq Mean Sq F value Pr(>F)
time 3 0.08205 0.027351 42.80 2.68e-14 ***
group:time 6 0.00272 0.000454 0.71 0.643
Residuals 54 0.03451 0.000639
If I run TukeyHSD:
TukeyHSD(anova)
Error in UseMethod("TukeyHSD") : no applicable method for 'TukeyHSD' applied to an object of class "c('aovlist', 'listof')"
Is anything wrong with my command or dataset? I know this is a very primary question... but if anyone could help me, it would be appreciated! Thanks.
Answer
What you're dealing with is repeated measures ANOVA, and you need to do the correct post-hoc testing for that. See below links for extra info:
Post hoc tests with ezANOVA output
Post hoc test after ANOVA with repeated measures using R
[R] Tukey HSD (or other post hoc tests) following repeated measures ANOVA
I think you're best off building a linear mixed-effects model with this specified error structure, as suggested in above links. Here is an artifical example dataset close to yours and post-hoc testing from multcomp-package for the model built using nlme-package:
set.seed(1)
dat <- cbind(expand.grid(time = paste("time", 1:3, sep=""), group = paste("task", 1:3, sep=""), subject = paste("S", 1:20, sep="")), response = rnorm(3*3*20))
# Add task1-specific effect (== task1.timeANY)
dat$response <- dat$response + as.numeric(dat$group=="task1")
# Extra effect in the last timepoint of task1 (== task1.time3)
dat$response <- dat$response + as.numeric(dat$group=="task1")*as.numeric(dat$time=="time3")
# Randomness specific for each subject
dat$response <- dat$response + rep(rnorm(20), each=3)
dat$grtim <- interaction(dat$group, dat$time)
# Interaction term specified above
#> head(dat)
# time group subject response grtim
#1 time1 task1 S1 -0.85777723 task1.time1
#2 time2 task1 S1 -0.04768010 task1.time2
#3 time3 task1 S1 -0.06695203 task1.time3
#4 time1 task2 S1 2.57917637 task2.time1
#5 time2 task2 S1 1.31340334 task2.time2
#6 time3 task2 S1 0.16342719 task2.time3
# Reason why TukeyHSD-function fails:
#anova = aov(response~(group*time)+Error(subject/time),dat)
#summary(anova)
#TukeyHSD(anova)
#Error in UseMethod("TukeyHSD") :
# no applicable method for 'TukeyHSD' applied to an object of class "c('aovlist', 'listof')"
#> class(anova)
#[1] "aovlist" "listof"
require(nlme)
# Below call does not work for glht, thus we created the interaction term in the data frame
#model <- lme(response ~ group*time, random = ~ 1 | subject / time, dat)
model <- lme(response ~ grtim, random = ~ 1 | subject / time, dat)
require(multcomp)
summary(glht(model, linfct=mcp(grtim="Tukey")), test = adjusted(type = "bonferroni"))
This outputs quite a long list of possible combinations, but we notice that task1, especially task1.time3, is quite different from the rest as expected:
Simultaneous Tests for General Linear Hypotheses
Multiple Comparisons of Means: Tukey Contrasts
Fit: lme.formula(fixed = response ~ grtim, data = dat, random = ~1 |
subject/time)
Linear Hypotheses:
Estimate Std. Error z value Pr(>|z|)
task2.time1 - task1.time1 == 0 -0.66574 0.40907 -1.627 1.000000
task3.time1 - task1.time1 == 0 -0.21758 0.40907 -0.532 1.000000
task1.time2 - task1.time1 == 0 0.46382 0.40907 1.134 1.000000
task2.time2 - task1.time1 == 0 -0.63987 0.40907 -1.564 1.000000
task3.time2 - task1.time1 == 0 -0.86698 0.40907 -2.119 1.000000
task1.time3 - task1.time1 == 0 1.17238 0.40907 2.866 0.149667
task2.time3 - task1.time1 == 0 -1.15241 0.40907 -2.817 0.174433
task3.time3 - task1.time1 == 0 -0.70811 0.40907 -1.731 1.000000
task3.time1 - task2.time1 == 0 0.44816 0.40907 1.096 1.000000
task1.time2 - task2.time1 == 0 1.12956 0.40907 2.761 0.207272
task2.time2 - task2.time1 == 0 0.02587 0.40907 0.063 1.000000
task3.time2 - task2.time1 == 0 -0.20124 0.40907 -0.492 1.000000
task1.time3 - task2.time1 == 0 1.83812 0.40907 4.493 0.000252 ***
task2.time3 - task2.time1 == 0 -0.48667 0.40907 -1.190 1.000000
task3.time3 - task2.time1 == 0 -0.04237 0.40907 -0.104 1.000000
task1.time2 - task3.time1 == 0 0.68140 0.40907 1.666 1.000000
task2.time2 - task3.time1 == 0 -0.42229 0.40907 -1.032 1.000000
task3.time2 - task3.time1 == 0 -0.64940 0.40907 -1.587 1.000000
task1.time3 - task3.time1 == 0 1.38996 0.40907 3.398 0.024451 *
task2.time3 - task3.time1 == 0 -0.93483 0.40907 -2.285 0.802723
task3.time3 - task3.time1 == 0 -0.49053 0.40907 -1.199 1.000000
task2.time2 - task1.time2 == 0 -1.10369 0.40907 -2.698 0.251098
task3.time2 - task1.time2 == 0 -1.33080 0.40907 -3.253 0.041077 *
task1.time3 - task1.time2 == 0 0.70856 0.40907 1.732 1.000000
task2.time3 - task1.time2 == 0 -1.61623 0.40907 -3.951 0.002802 **
task3.time3 - task1.time2 == 0 -1.17193 0.40907 -2.865 0.150188
task3.time2 - task2.time2 == 0 -0.22711 0.40907 -0.555 1.000000
task1.time3 - task2.time2 == 0 1.81225 0.40907 4.430 0.000339 ***
task2.time3 - task2.time2 == 0 -0.51254 0.40907 -1.253 1.000000
task3.time3 - task2.time2 == 0 -0.06824 0.40907 -0.167 1.000000
task1.time3 - task3.time2 == 0 2.03936 0.40907 4.985 2.23e-05 ***
task2.time3 - task3.time2 == 0 -0.28543 0.40907 -0.698 1.000000
task3.time3 - task3.time2 == 0 0.15887 0.40907 0.388 1.000000
task2.time3 - task1.time3 == 0 -2.32479 0.40907 -5.683 4.76e-07 ***
task3.time3 - task1.time3 == 0 -1.88049 0.40907 -4.597 0.000154 ***
task3.time3 - task2.time3 == 0 0.44430 0.40907 1.086 1.000000
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- bonferroni method)