如何使用Anova命令执行Tukey HSD测试(car包)

Question

我正在处理一个不平衡的设计/样本，并且最初学习了aov()。我现在知道，对于我的方差分析测试，我需要使用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测试？或者相反，使用Type III计算我的方差分析，并且仍然能够运行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函数括在圆括号中将导致它打印结果。请参见下面的示例。

一般来说，我建议使用flexible emmeans (née lsmeans)或multcomp包来满足您的所有后期比较需求。对于执行mean separations on interactions或examining contrasts among treatments，emmeans特别有用。编辑:注意:我是这些页面的作者。

对于不平衡的设计，您可能需要报告E.M. (或L.S.)均值而不是算术均值。参见SAEPER: What are least square means?。编辑:请注意，我是这个页面的作者。请注意，在下面的示例中，emmeans报告的边际平均值与HSD.test报告的边际平均值不同。

还要注意，glht中的" Tukey“与Tukey HSD或Tukey调整后的比较无关；它只是为所有成对测试设置对比度，如输出所示。

但是，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()的输出，必须严格使用MYaseen208答案中提供的逻辑来构造lm() (或者，如果模型使用aov()，也可以)：

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

希望这篇文章能帮助那些没有从HSD.test()得到正确输出的人。