Как сделать тест Tukey HSD с командой Anova (автомобильный пакет)

Я имею дело с несбалансированным дизайном/образцом и изначально узнал aov(). Теперь я знаю, что для тестов ANOVA мне нужно использовать Type of Squares, который включает использование фикчирования с использованием lm() вместо использования aov().

Проблема заключается в получении пост-hoc-тестов (в частности, Tukey HSD) с использованием lm(). Все проведенные исследования говорят, что использование simint в пакете multcomp будет работать, но теперь, когда он обновил эту команду, похоже, не будет доступным. Похоже, он полагается также на то, чтобы пройти aov() для вычисления.

По сути, все тесты Tukey HSD, которые я нашел для R, предполагают, что вы используете aov() для сравнения, а не для lm(). Чтобы получить сумму квадратов типа III, которые мне нужны для несбалансированного дизайна, я должен использовать:

mod<-lm(Snavg~StudentEthnicity*StudentGender)

Anova(mod, type="III")

Как использовать Tukey HSD-тест с модом, используя lm()? Или, наоборот, вычислить мой ANOVA с помощью Type III и все еще иметь возможность запускать тест Tukey HSD?

Спасибо!

Ответ 1

Попробуйте HSD.test в agricolae

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

Ответ 2

Я застрял в той же проблеме, когда HSD.test ничего не печатал. Вы должны поместить console=TRUE внутрь функции, чтобы она печаталась автоматически.

Например:

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

Ответ 3

Я нашел HSD.test() также очень тщательным образом о том, как вы построили модель lm() или aov(), которую вы используете для нее.

Не удалось получить данные из HSD.test() с моими данными, когда я использовал следующую идею кодирования для lm():

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

Ответ 4

В качестве первоначальной заметки, если она не была изменена, чтобы получить правильные результаты для суммы квадратов типа iii, вам необходимо установить контрастное кодирование фактор-переменных. Это можно сделать внутри вызова lm или с помощью options. В приведенном ниже примере используется options.

Я бы с осторожностью относился к использованию HSD.test и аналогичных функций с несбалансированными проектами, если только документация не обращается к их использованию в этих ситуациях. В документации для TukeyHSD упоминается, что она настраивается для "слегка несбалансированных" конструкций. Я не знаю, работает ли HSD.test по-разному. Вам нужно будет проверить дополнительную документацию для пакета или исходную ссылку, приведенную для этой функции.

В качестве побочного примечания, заключая в круглые скобки всю функцию HSD.test, он заставит ее распечатать результаты. См. Пример ниже.

В общем, я бы рекомендовал отказаться от Tukey.HSD и подобных функций и использовать гибкие пакеты emmeans (née lsmeans) или multcomp для всех ваших потребностей в постсоветском сопоставлении. emmeans особенно полезен для выполнения средних разделителей на взаимодействиях или для изучение контрастов между лечением.

С несимметричным дизайном вы можете сообщать о средствах E.M. (или L.S.) вместо арифметических средств. См. SAEPER: что такое наименьшее квадратное значение?. Обратите внимание, что в приведенном ниже примере маргинальные значения, сообщенные emmeans, отличаются от тех, которые указаны HSD.test.

Также обратите внимание, что "Tukey" в glht не имеет никакого отношения к сравнениям Tukey HSD или Tukey; он просто устанавливает контрасты для всех парных тестов, как говорит вывод.

Однако функции adjust="tukey" in emmeans означают использование сравнений, скорректированных Tukey, как говорит результат.

Следующий пример частично адаптирован из ARCHBS: односторонняя Anova.

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")}
if(!require(multcompView)){install.packages("multcompView")}
library(emmeans)
library(multcompView)

marginal = emmeans(model,
                   ~ cyl.f)

pairs(marginal, adjust="tukey")

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"))

multcomp::cld(mc)