r - R 中 TukeyHSD 中的错误

Question

我正在研究混合设计方差分析，并希望运行 TukeyHSD 进行事后测试。
我不断收到错误，“使用方法错误(“TukeyHSD”):
没有适用于 'TukeyHSD' 的方法应用于类 "c('aovlist', 'listof')"的对象。

我在一些 F&Q 网站上看到了有关此错误的几个问题，但我仍然找不到解决方案。

我的数据集如下所示:

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

然后，方差分析结果是:

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                     

如果我运行 TukeyHSD:

TukeyHSD(anova)
Error in UseMethod("TukeyHSD") : no applicable method for 'TukeyHSD' applied to an object of class "c('aovlist', 'listof')"

我的命令或数据集有什么问题吗？
我知道这是一个非常重要的问题……但如果有人能帮助我，将不胜感激!谢谢。

Answer 1

您正在处理的是重复测量方差分析，您需要为此进行正确的事后测试。有关其他信息，请参阅以下链接:

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

我认为您最好按照上述链接中的建议构建具有此指定错误结构的线性混合效应模型。这是一个接近你的人工示例数据集，以及来自 multcomp-package 的事后测试，用于使用 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"))

这会输出一长串可能的组合，但我们注意到 task1，尤其是 task1.time3，与预期的其他部分完全不同:

         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)