29 января 2013

Однофакторный дисперсионный анализ: введение



Рассмотренный ранее t-критерий Стьюдента (равно как и его непараметрические аналоги) предназначен для сравнения исключительно двух совокупностей. Однако часто он неверно используется для попарного сравнения большего количества групп (рис. 1), что вызывает т.н. эффект множественных сравнений (англ. multiple comparisons; Гланц 1999, с. 101-104). Об этом эффекте и о том, как с ним бороться, мы поговорим позднее. В этом же сообщении я опишу принципы  однофакторного дисперсионного анализа, как раз предназначенного для одновременного сравнения средних значений двух и более групп. Принципы дисперсионного анализа (англ. analysis of variance, ANOVA) были разработаны в 1920-х гг. сэром Рональдом Эйлмером Фишером (англ. Ronald Aylmer Fisher) - "гением, едва не в одиночку заложившим основы современной статистики" (Hald 1998).



Рис. 1. Пример неверного использования критерия Стьюдента для попарных сравнений трех групп - А, B и  C.

Может возникнуть вопрос: почему метод, используемый для сравнения средних значений, называется дисперсионным анализом? Все дело в том, что при установлении разницы между средними значениями мы в действительности сравниваем дисперсии анализируемых совокупностей. Однако обо всем по порядку...

Постановка задачи

Рассмотренный ниже пример заимствован из книги Maindonald & Braun (2010). Имеются данные о весе томатов (все растение целиком; weight, в кг), которые выращивали в течение 2 месяцев при трех разных экспериментальных условиях (trt, от  treatment)  - на воде (water), в среде с добавлением удобрения (nutrient), а также в среде с добавлением удобрения и гербицида 2,4-D (nutrient+24D):

# Создадим таблицу с данными:
tomato <-
  data.frame(weight=
               c(1.5, 1.9, 1.3, 1.5, 2.4, 1.5, # water
                 1.5, 1.2, 1.2, 2.1, 2.9, 1.6, # nutrient
                 1.9, 1.6, 0.8, 1.15, 0.9, 1.6), # nutrient+24D
             trt = rep(c("Water", "Nutrient", "Nutrient+24D"),
                       c(6, 6, 6)))
 
# Просмотрим результат:
weight
   weight          trt
1    1.50        Water
2    1.90        Water
3    1.30        Water
4    1.50        Water
5    2.40        Water
6    1.50        Water
7    1.50     Nutrient
8    1.20     Nutrient
9    1.20     Nutrient
10   2.10     Nutrient
11   2.90     Nutrient
12   1.60     Nutrient
13   1.90 Nutrient+24D
14   1.60 Nutrient+24D
15   0.80 Nutrient+24D
16   1.15 Nutrient+24D
17   0.90 Nutrient+24D
18   1.60 Nutrient+24D

Переменная trt представляет собой фактор с тремя уровнями. Для более наглядного сравнения экспериментальных условий в последующем, сделаем уровень "water" базовым (англ. reference), т.е. уровнем, с которым R будет сравнивать все остальные уровни. Это можно сделать при помощи функции relevel():

tomato$trt <- relevel(tomato$trt, ref = "Water")

Чтобы лучше понять свойства имеющихся данных, визуализируем их при помощи одномерной диаграммы рассеяния (рис. 2):

attach(tomato)
stripchart(weight ~ trt, xlab = "Вес, кг", ylab = "Условия")

Рис. 2. Результаты измерений веса растений томатов, выращенных при разных экспериментальных условиях.

Из рис. 2 видно, что измеренные значения веса растений достаточно близки для всех трех экспериментальных условий, хотя и есть некоторая тенденция к снижению веса в группе "Nutrient+24D". Это визуальное впечатление подтверждается также соответствующими групповыми средними:

tapply(weight, trt, mean)
       Water     Nutrient Nutrient+24D 
    1.683333     1.750000     1.325000

Подлежащую проверке нулевую гипотезу можно сформулировать так: исследованные условия выращивания растений не оказывают никакого влияния на вес последних. Другими словами, нулевая гипотеза утверждает, что наблюдаемые различия между групповыми средними несущественны и вызваны влиянием случайных факторов (т.е. в действительности все  полученные измерения веса растений происходят из одной нормально распределенной генеральной совокупности):

\[H_0: \mu_1 = \mu_2 = \mu_3\]

Подчеркнем еще раз, что рассматриваемый пример соответствует случаю однофакторного дисперсионного анализа: изучается действие одного фактора - условий выращивания (с тремя уровнями - Water, Nutrient и Nutrient+24D) на интересующую нас переменную-отклик - вес растений.

К сожалению, исследователь почти никогда не имеет возможности изучить всю генеральную совокупность. Как же нам тогда узнать, верна ли приведенная выше нулевая гипотеза, располагая только выборочными данными? Мы можем сформулировать этот вопрос иначе: какова вероятность  получить наблюдаемые различия между групповыми средними, извлекая случайные выборки из одной нормально распределенной генеральной совокупности? Для ответа на этот вопрос на нам потребуется статистический критерий, который количественно характеризовал  бы величину различий между сравниваемыми группами.




Перед тем, как сконструировать такой критерий, зададимся еще одним вопросом: что заставляет нас, взглянув на несколько графически представленных выборок, думать, что различия между ними неслучайны?

На рис. 3 к исходным данным добавлена еще одна группа - из точек, отражающих соответствующие выборочные средние (Means) .
Рис. 3.  То же, что рис. 2, но с добавлением точек, отражающих средние значения в каждой экспериментальной группе (Means).

Теперь (исключительно с целью продемонстрировать принцип!) несколько изменим исходные данные (рис. 3):
Рис. 4. То же, что рис. 3, но с искусственно измененными исходными данными.

Сравните рис. 4 с рис. 3 - что изменилось? Группы точек, отражающих экспериментальные данные, оказались значительно раздвинутыми вдоль оси X. Результатом этого стало также расхождение групповых средних (Means). Теперь, глядя на рис. 4, почти любой скажет, что экспериментальные группы различаются по весу растений. Почему? Сравните разброс значений внутри экспериментальных групп с разбросом трех групповых средних: разброс групповых средних на рис. 4 в целом превышает разброс значений в экспериментальных группах (тогда как на рис. 3 мы имели обратную ситуацию).

Следовательно, для оценки различий между группами следует каким-то образом сравнить разброс групповых средних с разбросом значений внутри групп. Это ключевая идея дисперсионного анализа, уяснив которую, вы не будете испытывать трудности с пониманием излагаемого ниже материала.


Две оценки дисперсии при дисперсионном анализе

Итак, чем больше разброс выборочных средних и чем меньше разброс значений внутри групп, тем меньше вероятность того, что наши группы являются случайными выборками из одной совокупности. Дисперсию генеральной совокупности можно оценить двумя способами. С одной стороны, оценкой дисперсии генеральной совокупностью будет дисперсия, вычисленная для каждой группы. Такая оценка не будет зависеть от различий групповых средних. С другой стороны, при верной нулевой гипотезе (см. выше) разброс  групповых средних тоже позволит оценить дисперсию генеральной совокупности. Очевидно, что такая оценка уже будет зависеть от различий между группами.

Если экспериментальные группы - это случайные выборки из одной и той же нормально распределенной генеральной совокупности, то оба способа оценки генеральной дисперсии должны давать примерно одинаковые результаты. Соответственно, если эти оценки действительно оказываются близки, то мы не можем отвергнуть нулевую гипотезу. И наоборот: если разница между этими оценками оказывается существенной, мы можем принять альтернативную гипотезу: маловероятно, что мы получили бы наблюдаемые различия между группами, если бы они были просто случайными выборками из одной нормально распределенной генеральной совокупности.

Перейдем к вычислениям. Пусть \(x_{ij}\) обозначает наблюдение \(j\) в группе \(i\) (например, \(x_{13}\) будет третьим наблюдением из первой группы), \(\bar{x_i}\) - среднее значение в группе \(i\), a \(\bar{X}\) - общее среднее значение (рассчитанное по всем имеющимся наблюдениям). Тогда каждое наблюдение мы можем разложить на следующие составляющие:

\[x_{ij} = \bar{X} + (\bar{x_i}- \bar{X}) + (x_{ij} - \bar{x_i})\]

где \((\bar{x_i}- \bar{X})\) - отклонения групповых средних от общего среднего значения, а \((x_{ij} - \bar{x_i})\) - отклонения отдельных наблюдений от среднего значения группы, к которой они принадлежат.

Тогда разброс наблюдений внутри групп можно рассчитать как

\[SS_W = \sum_{i}\sum_{j}(x_{ij} - \bar{x_i})^2\]

а разброс между группами (разброс групповых средних) как

\[SS_B = \sum_{i}\sum_{j}(\bar{x_i}- \bar{X})^2\]

В приведенных выражениях буквы W и B соответствуют английским словам "within" (внутри) и "between" (между). Нормализовав \(SS_W\) и \(SS_B\) по их соответствующим степеням свободы, получим внутри- и межгрупповую дисперсии, о которых мы говорили выше:

\[MS_W = SS_W/(N - k)\]
\[MS_B = SS_B/(k - 1)\]

В приведенных двух выражениях N - общее число наблюдений, а k - число сравниваемых групп. Аббревиатура MS означает "mean squares" ("средние квадраты"; имеются в виду усредненные суммы квадратов отклонений, \(SS_W\) и \(SS_B\)) и часто встречается в результатах дисперсионного анализа, выдаваемых статистическими программами.

Нормализация \(SS_W\) и \(SS_B\) по числу степеней свободы позволяет получить сравнимые величины. Формально внутри- и межгрупповые дисперсии сравниваются при помощи F-критерия, или критерия Фишера (обратите внимание: в числителе всегда находится межгрупповая дисперсия, \(MS_B\)):

\[F = MS_B/MS_W\]

Очевидно, что чем ближе F к 1, тем меньше у нас оснований утверждать, что внутри- и межгрупповая дисперсии различаются. Иными словами, у нас нет оснований отклонить сформулированную выше нулевую гипотезу. Если же F значительно выше 1, нулевую гипотезу можно отклонить. Возникает вопрос: начиная с какой именно величины F нулевую гипотезу можно отвергать?


Критическое значение F-критерия

Критическое значение F-критерия определяется желаемым уровнем значимости и свойствами F-распределения, форма которого полностью задается меж- и внутригрупповым степенями свободы. Так, для нашего примера, межгрупповое число степеней свободы составляет \(df_B = k - 1 =  3 - 1 = 2\), а внутригрупповое - \(df_W = N - k = 18 - 3 = 15\). Внешний вид F-распределения при этих \(df_B\) и \(df_W\) представлен на рис. 5. Вертикальная линия на этом рисунке соответствует 3.682 - критическому значению F при \(\alpha = 0.05\). Если F-значение, рассчитанное по экспериментальным данным, превышает критическое значение, мы можем отклонить нулевую гипотезу об отсутствии эффекта изучаемого фактора.


Выполнение дисперсионного анализа в R

Дисперсионный анализ в R можно выполнить при помощи базовых функций aov() и lm(). В этом сообщении мы рассмотрим только функцию aov(). Для нашего примера получаем:

summary(aov(weight ~ trt, data = tomato))
            Df Sum Sq Mean Sq F value Pr(>F)
trt          2  0.627  0.3135   1.202  0.328
Residuals   15  3.912  0.2608

В приведенных результатах строка, обозначенная как trt, соответствует источнику дисперсии в данных, связанному с действием изучаемого экспериментального фактора - условий выращивания растений. Строка, обозначенная как Residuals, характеризует внутригрупповую дисперсию (ее еще называют шумовой или остаточной дисперсией - в том смысле, что она не может быть объяснена влиянием экспериментального фактора). Столбец Sum Sq содержит \(SS_B\) и \(SS_W\), а столбец Mean Sq -  меж- и внутригрупповую дисперсии, \(MS_B\) и \(MS_W\). В столбце F value представлено рассчитанное по имеющимся данным значение F-критерия. Наконец, в столбце Pr(>F) представлена вероятность получить F-значение, равное или превышающее то значение, которое мы в действительности рассчитали по имеющимся выборочным данным (при условии, что нулевая гипотеза верна). Как видим, эта вероятность достаточно высока. Во всяком случае, она превышает 5%-ный уровень значимости, в связи с чем мы заключаем, что нулевая гипотеза верна. Таким образом, с достаточно высокой степенью уверенности мы можем утверждать, что экспериментальные условия не оказали существенного влияния на вес растений.

Как следует из названия, данное сообщение является введением в дисперсионный анализ и его выполнение при помощи R. В последующих сообщениях будут обсуждены такие вопросы, как условия применимости параметрического дисперсионного анализа и способы их проверки, множественные сравнения групп, расчет необходимого числа наблюдений для дисперсионного анализа, и др. "Не переключайтесь!"


Литература:

  • Гланц С (1999) Медико-биологическая статистика. - М.: Практика
  • Hald A (1998) A History of Mathematical Statistics. - New York: Wiley
  • Maindonald JH, Braun WJ (2010) Data Analysis and Graphics Using R. - Cambridge University Press


12 комментариев :

valinn комментирует...

Очень понятно объясняете. Получил удовольствие от чтения. Большое спасибо!

Анонимный комментирует...

спасибо огромнейшее!
мечты-мечты, но если бы Вы как-нибудь вдруг нашли время сделать подобное по ковариационному анализу. вот типа такого, но на русском http://r-eco-evo.blogspot.ru/2011/08/comparing-two-regression-slopes-by.html , было бы просто супер.

Сергей Мастицкий комментирует...

Всем спасибо за комментарии!
ANCOVA в планах, причем в ближайших. Так что подписываейтесь на рассылку, чтобы не пропустить :)

Анонимный комментирует...

Ваш блог - "настольная книга", я тут постоянно "пасусь" (да и подписан на Вашу R группу вконтакте), так что не пропущу. Еще раз спасибо!

Анонимный комментирует...

Простите, что беспокою. Вы не могли бы уточнить (чуть расширить или перефразировать) вот этот абзац:"Дисперсию генеральной совокупности можно оценить двумя способами. С одной стороны, оценкой дисперсии генеральной совокупностью будет дисперсия, вычисленная для каждой группы. Такая оценка не будет зависеть от различий групповых средних. С другой стороны, при верной нулевой гипотезе (см. выше) разброс групповых средних тоже позволит оценить дисперсию генеральной совокупности. Очевидно, что такая оценка уже будет зависеть от различий между группами."
все-таки не очень понятно "..будет дисперсия, вычисленная для каждой группы". Если мы исходим из нулевой гипотезы (наши выборки из одной ген совокупности), тогда, наверное, "..будет дисперсия, вычисленная для всех групп". Нет? возможно, я что-то недопониимаю. Спасибо заранее за Ваше терпение!

Сергей Мастицкий комментирует...

"..будет дисперсия, вычисленная для каждой группы" - речь идет о внутригрупповой дисперсии. Для каждой группы считается своя дисперсия, после чего эти групповые дисперсии усредняются (рассмотрите еще раз хорошенько формулы для SSw и MSw)

В свою очередь межгрупповая дисперсия рассчитывается по групповым средним (относительно общего среднего значения, рассчитанного по все имеющимся наблюдениям) - см. формулы для SSB и MSB.

Анонимный комментирует...

Спасибо!

zeleniy комментирует...

Сергей, расскажите пожалуйста, как вы дорисовали к stripchart'у средние значения means для каждой выборки соответствущими цветами. Не в фотошопе, надеюсь :)

Sergey Mastitsky комментирует...

Нет, до фотошопа дело не дошло :)
Функция striphchart(), к сожалению, действительно не позволяет избирательно залить разным цветом точки, принадлежащие одной группе (в нашем случае Means). Поэтому пришлось ее немного обмануть - сначала отобразить точки так, как эта функция делает автоматически, а потом при помощи функции points() поверх нужных точек нарисовать новые с желаемыми цветами. Весь код выглядел примерно так:

# Таблица со средними значениями:
Means <- data.frame(weight = as.numeric(tapply(tomato$weight,
tomato$trt, mean)),
trt = rep("Means", 3))

# Добавляем таблицу Means к таблице tomato:
tomato <- rbind(tomato, Means)

# Изменяем базовый уровень фактора trt на Water:
tomato$trt <- relevel(tomato$trt, ref = "Water")

# Рисуем исходную диаграмму (все точки из группы Means будут при этом
# автомтически залиты черным цветом):
stripchart(weight ~ trt, data = tomato, pch = 19,
col = c("blue", "red", "black"))

# Добавляем поверх точек из группы Means точки нужного цвета
# (эти новые точки имеют координаты (средняя_1, 4), (средняя_2, 4)
# и (средняя_3, 4); их цвет прописываем "вручную"):
points(x = Means$weight, y = c(4, 4, 4), pch = 19,
col = c("red", "black", "blue"))

zeleniy комментирует...

спасибо :)

Andrey ILIN комментирует...

Здорово!

Анонимный комментирует...

спасибо огромное !!!!!!!

Отправить комментарий