29 августа 2013

Классические методы статистики: дисперсионный анализ по Краскелу-Уоллису



Как было отмечено ранее, важными условиями применимости классического однофакторного дисперсионного анализа являются нормальность распределения зависимой переменной и однородность (гомоскедастичность) дисперсий во всех сравниваемых группах. В случаях, когда наблюдается существенное нарушение этих условий и ситуацию не получается исправить путем трансформации исходных значений анализируемой переменной (см. Box & Cox 1964), решением может стать применение дисперсионного анализа по Краскелу-Уоллису (англ. Kruskal-Wallis ANOVA by ranks или Kruskal-Wallis rank sum test; см. также оригинальную статью c описанием метода: Kruskal & Wallis 1952). В русскоязычной литературе для этого метода используются также названия "критерий Крускала-Уоллиса", "Н-критерий Крускала-Уоллиса" и даже иногда "критерий Крускала-Валлиса". В этом сообщении я покажу, как тест Краскела-Уоллиса выполняется в программе R.


Немного теории

Дисперсионный анализ по Краскелу-Уоллису относится к группе непараметрических методов статистики. Это значит, что при выполнении соответствующих расчетов параметры того или иного вероятностного распределения (например, нормального) никак не задействованы. Вместо этого используются ранги исходных значений и их суммы в сравниваемых группах. В частности, метод Краскела-Уоллиса основан на вычислении т.н. H-критерия:

\[ H=\frac{12}{N(N+1)}\sum_{i=1}^{m}\frac{R_i^{2}}{n_i}-3(N+1),\]

где \(n_i\) - число наблюдений в группе \(i\), \(N = \sum_{i=1}^{m}n_i \) - общее число наблюдений во всех \(m\) группах, а \(R_i\) - сумма рангов наблюдений в группе \(i\). Ранг представляет собой порядковый номер конкретного наблюдения в ряду упорядоченных по возрастанию (или убыванию - не важно) наблюдений (см. также статью по критерию Уилкоксона). Чем больше значение Н-критерия, тем больше у нас оснований отклонить нулевую гипотезу об отсутствии разницы между сравниваемыми группами. Если рассчитанное по выборочным данным значение Н превышает определенное критическое значение, нулевая гипотеза отклоняется. Критическое значение определяется с учетом принятого уровня значимости и числа степеней свободы; в частности, при \(m > 5\) H-критерий сравнивается с критическими значениями критерия хи-квадрат для числа степеней свободы \(m - 1\). При меньшем числе сравниваемых групп вносятся определенные поправки (например, поправка Имана-Давенпорта).

Интересно, что если бы мы выполнили обычный дисперсионный анализ на основе ранговых номеров исходных значений анализируемой переменной, то результат совпал бы результатом теста Краскела-Уоллиса (Zar 1999). Отсюда использование "дисперсионного анализа" в названии метода Краскела-Уоллиса (см. выше). Кроме того, при наличии двух сравниваемых групп, тест Краскела-Уоллиса будет идентичен тесту Манна-Уитни.

Если бы анализируемые данные удовлетворяли условиям нормальности и однородности групповых дисперсий, то статистическая мощность теста Краскела-Уоллиса в отношении таких данных составила бы примерно 95% от обычного параметрического дисперсионного анализа. Однако при нарушении этих условий мощность тест Краскела-Уоллиса может оказаться даже выше, чем у обычного дисперсионного анализа (Zar 1999).

В теории, для расчета несмещенных оценок Н-критерия Краскела-Уоллиса значения анализируемой переменной должны иметь одинаковый разброс и форму распределения во всех сравниваемых группах. Однако на практике нарушение этих условий мало сказывается на качестве выводов, получаемых при помощи Н-критерия, и ими обычно пренебрегают (Zar 1999).


Функция kruskal.test()


В R дисперсионный анализ по Краскелу-Уоллису выполняется при помощи функции kruskal.test(). В качестве примера (уже в который раз!) используем данные эксперимента по изучению эффективности 6 разных инсектицидов (таблица InsectSprays с этими данными входит в состав базовой версии R; см. также здесьздесь и здесь). Каждым из этих средств обработали по 12 растений, после чего подсчитали количество выживших насекомых. Нулевая гипотеза, которую мы проверим при помощи теста Краскела-Уоллиса, заключается в том, что исследованные инсектициды не различаются по эффективности (т.е. наблюдаемые различия в групповых медианных значениях числа выживших насекомых совершенно случайны).


kruskal.test(count ~ spray, data = InsectSprays)
 
 Kruskal-Wallis rank sum test
 
data:  count by spray
Kruskal-Wallis chi-squared = 54.6913, df = 5, p-value = 1.511e-10

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

Следует подчеркнуть, что подобно классическому дисперсионному анализу, тест Краскела-Уоллиса позволяет сделать заключение только следующего вида: либо "сравниваемые группы статистически значимо различаются" (например, при Р < 0.05), либо "статистически значимых различий между группами нет" (например, при Р > 0.05). Ни один из этих методов сам по себе не позволяет сказать, где именно лежат различия. Чтобы выяснить это необходимо выполнить соответствующие апостериорные тесты (англ. post hoc tests). Этой теме будут посвящены следующие несколько сообщений.


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

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

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

Правильно я понимаю, "одинаковый разброс" - это и есть условие равенства групповых дисперсий?

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

Да, поскольку разброс (=вариацию, изменчивость) выражают при помощи дисперсии, можно говорить об условии однородности групповых дисперсий. Но, как отмечено в статье, невыполнение этого условия почти не оказывает влияния на конечный результат, поскольку метод основан на рангах, а ранги существенно "сглаживают" различия между отдельными наблюдениями. Например, сумма рангов в ряду (10, 20, 30) составила бы 1+2+3 = 6. Идентичную сумму рангов получим и для ряда (1, 100, 1000), где разброс (дисперсия) гораздо выше.

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

Сергей, ещё такой вопрос.
А что если мы имеем дело с двумя факторами, можно ли в этом случае применять метод Краскела-Уоллиса, как аналог многофакторного дисперсионного анализа? Или для двух факторного экспериментального плана существуют свои непараметрические аналоги?

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

Метод Краскела-Уоллиса предназначен только для ситуаций с одним фактором. При наличии двух и более факторов можно применить следующие подходы:

1) Преобразовать исходные данные в ранги (см. хелп-файл по функции rank) и потом выполнить обычный двухфакторный дисперсионный анализ по этим преобразованным данным.

2) Использовать перестановочные тесты. Подробнее см. книгу Шитикова и Розенберга (2013): http://adf.ly/TC4za

3) Тест Краскела-Уоллиса можно рассматривать как частный случай т.н. "proportional odds regression" (регрессионной модели пропорциональных рисков; поищите в Google - на английском языке по этой теме много материала). В этом случае исходную переменную преобразовывают в фактор с упорядоченными уровнями. При небольшом разбросе в данных, метками уровней могут служить исходные значения (т.е. уровней будет столько же, сколько уникальных значений исходной переменной). При большом разбросе, возможно, придется объединить несколько соседних значений в более крупные группы (что, конечно, приведет к утрате определенной части заключенной в данных информации).
Такой анализ в R можно выполнить при помощи нескольких функций:
- lrm() из пакета rms
- clm() из пакета ordinal
- polr() из пакета MASS (входит в базовый набор пакетов R)
Для грамотного использования перечисленных функций и правильной интерпретации получаемых результатов нужно будет почитать соответствующую литературу. В отношении анализа экологических данных могу посоветовать эту статью: http://adf.ly/UxIdw

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

"1) Преобразовать исходные данные в ранги (см. хелп-файл по функции rank) и потом выполнить обычный двухфакторный дисперсионный анализ по этим преобразованным данным."

Сергей,
имеет ли в таком случае значение, какую из функций применять - aov() или lm() - будет ли ранжирование данных как-то влиять на результат?
На сколько можно судить из Ваший ранних сообщений, в случае с несбалансированными выборками, дисперсионный анализ нужно рассматривать как линейную модель, т.е. применять функцию lm().

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

Логика остается прежней: если число наблюдений в группах примерно одинаково, можно использовать aov(), если нет - lm(). В целом, использование lm() будет более "безопасным".

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

Скажите, если у нас группирующая переменная - дихотомическая, а значения в группах измеряются по порядковой шкале (1-7), какой метод для поиска наличия/отсутствия связи лучше использовать?

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