如何使用 Anova 命令进行 Tukey HSD 测试（汽车包）

Question

我正在处理一个不平衡的设计/样品，最初学习了aov()。我现在知道，对于我的 ANOVA 测试，我需要使用 III 型平方和，这涉及使用 lm() 而不是使用 aov() 进行拟合。

问题在于使用lm() 进行事后测试（特别是 Tukey 的 HSD）。我所做的所有研究都表明在multcomp 包中使用simint 会起作用，但现在它已更新，该命令似乎不可用。它似乎也依赖通过aov()来计算。

基本上我为 R 找到的所有 Tukey HSD 测试都假设您使用 aov() 进行比较而不是 lm()。要获得我必须使用的不平衡设计所需的 III 型平方和：

mod<-lm(Snavg~StudentEthnicity*StudentGender)

Anova(mod, type="III")

如何使用 lm() 对我的 mod 进行 Tukey HSD 测试？或者相反，使用类型 III 计算我的 ANOVA 并且仍然能够运行 Tukey HSD 测试？

谢谢！

Answer 1

在agricolae 中尝试HSD.test

library(agricolae)
data(sweetpotato)
model<-lm(yield~virus, data=sweetpotato)
comparison <- HSD.test(model,"virus", group=TRUE,
main="Yield of sweetpotato\nDealt with different virus")

输出

Study: Yield of sweetpotato
Dealt with different virus

HSD Test for yield 

Mean Square Error:  22.48917 

virus,  means

      yield  std.err replication
cc 24.40000 2.084067           3
fc 12.86667 1.246774           3
ff 36.33333 4.233727           3
oo 36.90000 2.482606           3

alpha: 0.05 ; Df Error: 8 
Critical Value of Studentized Range: 4.52881 

Honestly Significant Difference: 12.39967 

Means with the same letter are not significantly different.

Groups, Treatments and means
a        oo      36.9 
ab       ff      36.33333 
 bc      cc      24.4 
  c      fc      12.86667 

Answer 2

作为初始说明，除非已更改，否则要获得类型 iii 平方和的正确结果，您需要为因子变量设置对比度编码。这可以在lm 调用或options 中完成。下面的例子使用options。

我会谨慎使用 HSD.test 和具有不平衡设计的类似功能，除非文档说明了它们在这些情况下的使用。 TukeyHSD 的文档提到它针对“轻度不平衡”设计进行了调整。我不知道HSD.test 的处理方式是否不同。您必须检查包的其他文档或函数引用的原始参考。

附带说明，将整个 HSD.test 函数括在括号中将导致它打印结果。请参见下面的示例。

一般来说，我建议您使用灵活的emmeans（néelsmeans）或multcomp 包来满足您的所有事后比较需求。 emmeans 对于执行mean separations on interactions 或examining contrasts among treatments 特别有用。 [编辑：请注意，我是这些页面的作者。]

对于不平衡设计，您可能希望报告 E.M.（或 L.S.）均值而不是算术均值。见SAEPER: What are least square means?。 [编辑：请注意，我是本页的作者。] 请注意，在下面的示例中，emmeans 报告的边际平均值与 HSD.test 报告的边际平均值不同。

还要注意glht 中的“Tukey”与 Tukey HSD 或 Tukey-adjusted 比较无关；正如输出所示，它只是为所有成对测试设置了对比。

但是，正如输出所示，emmeans 函数中的 adjust="tukey" 确实意味着使用 Tukey 调整的比较。

以下示例部分改编自ARCHBS: One-way Anova。

### EDIT: Some code changed to reflect changes to some functions
###  in the emmeans package

if(!require(car)){install.packages("car")}
library(car)
data(mtcars)
mtcars$cyl.f = factor(mtcars$cyl)
mtcars$carb.f = factor(mtcars$carb)

options(contrasts = c("contr.sum", "contr.poly"))

model = lm(mpg ~ cyl.f + carb.f, data=mtcars)

library(car)
Anova(model, type="III")

if(!require(agricolae)){install.packages("agricolae")}
library(agricolae)
(HSD.test(model, "cyl")$groups)

if(!require(emmeans)){install.packages("emmeans")}
library(emmeans)

marginal = emmeans(model,
                   ~ cyl.f)

pairs(marginal, adjust="tukey")

if(!require(multcomp)){install.packages("multcomp")}
library(multcomp)

cld(marginal, adjust="tukey", Letters=letters)


if(!require(multcomp)){install.packages("multcomp")}
library(multcomp)

mc = glht(model,
          mcp(cyl.f = "Tukey"))

summary(mc, test=adjusted("single-step"))

cld(mc)

Answer 3

我发现HSD.test() 对您构建用于它的lm() 或aov() 模型的方式也非常细致。

当我对lm() 使用以下编码思想时，我的数据没有来自HSD.test() 的输出：

    model<-lm(sweetpotato$yield ~ sweetpotato$virus)  
    out <- HSD.test(model,"virus", group=TRUE, console=TRUE)

输出只有：

    Name:  virus 
    sweetpotato$virus 

当对aov() 使用相同的逻辑时，输出同样糟糕

    model<-aov(sweetpotato$yield ~ sweetpotato$virus)

获取HSD.test() 的输出lm() （或者如果模型使用aov()） 必须使用 MYaseen208 答案中提供的逻辑严格构建：

    model <- lm(yield~virus, data=sweetpotato)

希望这可以帮助那些没有从 HSD.test() 获得正确输出的人。

Answer 4

我遇到了同样的问题，即 HSD.test 什么也没打印出来。您需要将console=TRUE 放入函数中，以便它自动打印出来。

例如：

HSD.test(alturacrit.anova, "fator", console=TRUE).
Hope it helps!