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)
