R: TUKEY после одностороннего ANOVA
Я провел односторонний ANOVA в R, но продолжаю получать сообщения об ошибках, когда я пытаюсь сделать Tukey post-hoc, чтобы увидеть, какие обработки отличаются друг от друга. (Я бы хотел, чтобы результаты оценивались (a, ab, b, bcd... и т. Д.)
Детали ДАННЫХ:
data = "abh2"
х = 6 процедур: "лечение"
y = показания влажности "влажный" (n=63 на обработку, всего =378)
Я запустил одностороннюю ANOVA:
anov <- anova(lm(moist~treatment, data=abh2))
. # РЕЗУЛЬТАТЫ указывают на то, что я могу перейти в специальное состояние (р<0,05):
Analysis of Variance Table
Response: moist
Df Sum Sq Mean Sq F value Pr(>F)
treatment 5 1706.3 341.27 25.911 < 2.2e-16 ***
Я выбрал Tukey HSD и попытался запустить его двумя способами, но получил сообщения об ошибках для обоих:
Встроенная функция R:
TukeyHSD(anov)
# ERROR : no applicable method for 'TukeyHSD' applied to an object of class "c('anova', 'data.frame')"
Агрикола пакет:
HSD.test(anov, "treatment", group=TRUE, console=TRUE)
# ERROR : Error in HSD.test(anov, "treatment", group = TRUE, console = TRUE) :
argument "MSerror" is missing, with no default
Я обнаружил, что MSerror был
1) "# Старая версия HSD.test()" (но я только что обновил пакет Agricolae)
2) MSerror<-deviance (модель) / df
Итак, я попробовал:
HSD.test(anov, "treatment", MSerror=deviance(moist)/5, group=TRUE, console=TRUE)
*but still* # ERROR: $ operator is invalid for atomic vectors
Может ли кто-нибудь помочь мне продвинуться отсюда? Это кажется довольно простой проблемой, но я потратил на это часы!
Большое спасибо:)
2 ответа
Спасибо за обратную связь с Энни-Клод, я понял, что R не распознает данные должным образом.
Я решил проблему с помощью этого кода, хотя:
library(agricolae)
model<-aov(moist~treatment, data=abh2)
out <- HSD.test(model,"treatment", group=TRUE,console=TRUE)
Похоже, что команда ANOVA, которую я первоначально использовал (из базового пакета R), просто не была понята тестом Тьюки, который я пытался выполнить впоследствии (с пакетом Agricolae).
Возьмите с собой сообщение для дома: попробуйте провести соответствующую серию анализов в одном пакете!
ps Для получения значений p:
summary(model)
они могут быть преобразованы в фрейм данных следующим образом:
as.data.frame (сводка (модель)[[1]])
Попробуйте указать ваше лечение как фактор с помощью следующего кода:
abh2$treatment <- factor(abh2$treatment)