Дисперсионный анализ в R: aov(), таблица ANOVA и критерий Тьюки

Почему дисперсионный анализ, а не гора t-критериев

t-критерий сравнивает два средних. При трёх и более группах возникает соблазн проверить t-критерием каждую пару — и это ровно ловушка. Каждый критерий на уровне 0.05 несёт 5-процентный риск ложноположительного результата, и риски накапливаются: при 4 группах это 6 попарных критериев и примерно 26-процентный шанс хотя бы одного ложного «значимого» результата; при 5 группах (10 критериев) — около 40%. Вы бы фабриковали открытия из чистого шума.

Дисперсионный анализ (ANOVA, ANalysis Of VAriance) исправляет это, задавая один вопрос с одним p-значением: все ли групповые средние одинаковы или хотя бы одно отличается? Несмотря на название, он сравнивает средние — просто делает это через анализ дисперсии: если групповые средние разбросаны сильнее, чем может объяснить шум внутри групп, происходит что-то настоящее.

Однофакторный дисперсионный анализ через aov()

PlantGrowth создан для этого: вес растений в контроле и при двух условиях обработки. Подгоните через aov(), а затем — это важно — напечатайте таблицу через summary():

fit <- aov(weight ~ group, data = PlantGrowth)
summary(fit)

Читайте формулу как «зависит ли weight от group?» Группирующая переменная должна быть фактором — PlantGrowth$group им уже является, но если ваши группы закодированы числами (дозы, номера партий), оберните их: aov(y ~ factor(dose), ...). Иначе aov() молча подгонит линию регрессии по кодам групп вместо сравнения групповых средних — неверная модель без всякого сообщения об ошибке.

Чтение таблицы дисперсионного анализа

Таблица состоит из двух строк — фактор и остатки — и пяти столбцов. Ячейка за ячейкой:

            Df Sum Sq Mean Sq F value Pr(>F)
group        2  3.766  1.8832   4.846 0.0159
Residuals   27 10.492  0.3886

• Df — число степеней свободы. Строка фактора получает число групп − 1 (3 группы → 2); строка остатков получает число наблюдений − число групп (30 − 3 = 27). Быстрая проверка того, что R увидел задуманный вами дизайн.
• Sum Sq — изменчивость, разделённая на две кучи. Строка group — межгрупповая куча: насколько групповые средние отстоят от общего среднего. Residuals — внутригрупповая куча: насколько растения варьируют вокруг среднего своей группы. Вместе они складываются в полную изменчивость данных.
• Mean Sq — каждая Sum Sq, делённая на своё Df, что превращает кучи изменчивости в сопоставимые величины на степень свободы. Остаточная Mean Sq (0.389) — это уровень шума.
• F value — отношение: Mean Sq строки group к Mean Sq строки Residuals (1.8832 / 0.3886 ≈ 4.85). Если бы все групповые средние были действительно равны, это отношение колебалось бы около 1. Чем больше F, тем сильнее межгрупповые различия обгоняют то, что может объяснить шум.
• Pr(>F) — p-значение: вероятность получить настолько большое F, если бы все три истинных средних были одинаковы. Здесь 0.016 — достаточно мало при обычном уровне 0.05, чтобы отвергнуть «все равны». Как всегда, это не вероятность истинности нулевой гипотезы, и оно ничего не говорит о том, какие группы различаются и насколько.

Последний пункт и есть ключевое ограничение: значимое F говорит «где-то есть какое-то различие». И ничего больше.

Какие группы различаются? TukeyHSD()

Следующему вопросу нужен апостериорный критерий. Критерий Тьюки (Honest Significant Difference) проверяет каждую пару, удерживая групповую частоту ошибок на уровне 5% сразу по всем сравнениям:

fit <- aov(weight ~ group, data = PlantGrowth)
TukeyHSD(fit)

По строке на пару, по четыре числа в строке:

• diff — оценённая разность средних (вторая названная группа минус первая).
• lwr, upr — 95-процентный групповой доверительный интервал для этой разности.
• p adj — p-значение, уже скорректированное на три сравнения.

Для PlantGrowth: trt2-trt1 показывает разность около 0.87 с p adj ≈ 0.012 и интервалом, не задевающим ноль, — обработка 2 обгоняет обработку 1. У обеих строк «обработка против контроля» (trt1-ctrl, trt2-ctrl) интервалы накрывают ноль, а p adj заметно выше 0.05 — ни одна обработка на этой выборке неотличима от контроля. Правило большого пальца зеркалит доверительные интервалы повсюду: интервал не содержит нуля ⇔ p adj ниже 0.05.

Обратите внимание, как общее F сказало «различие существует», а Тьюки локализует его ровно в одной паре, — двухшаговая структура и есть весь замысел: один честный общий критерий, а затем должным образом скорректированная попарная детективная работа.

Двухфакторный анализ: два фактора и их взаимодействие

При двух группирующих переменных один вызов проверяет обе, а также то, взаимодействуют ли они. ToothGrowth скрещивает тип добавки (supp) с дозой (0.5, 1, 2 мг — числовая, поэтому ей нужен factor()):

fit2 <- aov(len ~ supp * factor(dose), data = ToothGrowth)
summary(fit2)

supp * factor(dose) разворачивается в три эффекта, по строке таблицы на каждый:

• supp — в среднем по дозам, важен ли тип добавки? (Да: p ≈ 0.0002.)
• factor(dose) — в среднем по добавкам, важна ли доза? (Категорически да: p ничтожно мало.)
• supp:factor(dose) — взаимодействие: меняется ли эффект добавки в зависимости от дозы? Здесь p ≈ 0.022 — меняется. Апельсиновый сок обыгрывает аскорбиновую кислоту при низких дозах, но при 2 мг разрыв смыкается.

Значимое взаимодействие — предупреждающая наклейка на главных эффектах: утверждение «тип добавки важен» верно лишь в среднем, а среднее прячет зависящую от дозы историю. Когда взаимодействие значимо, описывайте комбинации (через TukeyHSD(fit2, "supp:factor(dose)") или таблицу групповых средних из групповых сводок), а не сообщайте одни только главные эффекты. Используйте + вместо * только когда вы намеренно хотите модель без взаимодействия.

Предпосылки и ранговый запасной вариант

Математика дисперсионного анализа опирается на три предпосылки, в порядке убывания договороспособности:

• Независимость — наблюдения не влияют друг на друга. Нарушение ничем не исправить задним числом; это свойство дизайна исследования.
• Равенство дисперсий по группам — остаточная Mean Sq является одной объединённой оценкой шума, поэтому группы должны быть примерно одинаково шумными. Однострочная проверка: bartlett.test(weight ~ group, data = PlantGrowth) (большое p-значение означает отсутствие свидетельств неравенства дисперсий).
• Примерно нормальные остатки — имеет наименьшее значение при сбалансированных планах и приличных размерах групп; окиньте взглядом гистограмму или ящик с усами остатков через residuals(fit).

Когда данные сильно скошены, порядковые или изъедены выбросами, ранговой альтернативой однофакторному дисперсионному анализу служит kruskal.test(weight ~ group, data = PlantGrowth) — многогрупповой родственник wilcox.test(). Это одна строка, меняющая часть мощности на устойчивость.

Что вы уносите с собой

• Дисперсионный анализ задаёт один вопрос о средних 3+ групп с одним p-значением — это лекарство от множественности t-критериев.
• fit <- aov(y ~ group, data = df); summary(fit) — и группирующая переменная должна быть фактором.
• Значение F — это межгрупповой сигнал, делённый на внутригрупповой шум; малое Pr(>F) означает «где-то есть какое-то различие», и ничего больше.
• TukeyHSD(fit) находит, какие пары различаются, со встроенной поправкой на множественные сравнения.
• a * b подгоняет два фактора плюс их взаимодействие; значимое взаимодействие означает, что главные эффекты сами по себе не рассказывают всю историю.
• Проверяйте равенство дисперсий (bartlett.test) и нормальность остатков; kruskal.test() — ранговый запасной вариант.

Дальше: от сравнения групповых средних к моделированию связи — линейная регрессия через lm().