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

Question

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

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

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

mod<-lm(Snavg~StudentEthnicity*StudentGender)

Anova(mod, type="III")

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

Спасибо!

Answer 1

Попробуйте HSD.test in 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 

Answer 2

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

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

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

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

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

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

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

Следующий пример частично адаптирован из 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)

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

Answer 4

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

Например:

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