20.05.2012

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

В одном из предыдущих сообщений я описал расчет критерия Стьюдента при помощи соответствующей функции, входящей в базовую комплектацию R - t.test(). Одно из важных условий корректного применения критерия Стьюдента состоит в том, что анализируемые выборки должны происходить из нормально распределенных генеральных совокупностей. В случаях, когда это условие не выполняется, вместо критерия Стьюдента следует использовать  его непараметрический аналог - критерий Уилкоксона (Wilcoxon rank test). Здесь следует сразу пояснить, что создатели системы R под названием "критерий Уилкоксона" (или "тест Уилкоксона") объединяют как собственно метод, предложенный Фрэнком Уилкоксоном (Frank Wilcoxon) в 1945 г., так и опубликованный несколько позднее (1947 г.) метод Манна-Уитни. Первый из этих методов обычно используется для сравнения двух связных выборок, тогда как второй предназначен для сравнения двух независимых выборок. Ниже я не буду разграничивать эти методы, придерживаясь подхода, принятого в системе R.


Одновыборочный критерий Уилкоксона

Этот вариант критерия (Wilcoxon signed rank test) служит для проверки нулевой гипотезы о том, что анализируемая выборка происходит из симметрично распределенной генеральной совокупности с центром в точке µ0. Алгоритм метода заключается в следующем:
  • µ0 отнимают от каждого выборочного значения;
  • получившиеся величины ранжируют по возрастанию, игнорируя знак;
  • ранговые номера со знаком + суммируют, получая величину V;
  • критерий V сравнивают с критическим значением для заданного уровня значимости и числа степеней свободы; альтернативный подход на этом шаге заключается в нахождении вероятности случайным образом получить значение V, равное или превышающее наблюдаемое значение, при условии истинности нулевой гипотезы.
Ниже при описании теста Уилкоксона я буду использовать примеры, рассмотреные ранее для иллюстрации t-теста Стьюдента. Это позволит нам сравнить результаты, получаемые при помощи обоих методов. В частности, обратимся к данным о суточном потреблении энергии (кДж/сутки) у 11 женщин (пример заимствован из книги Altman D. G. (1981) Practical Statistics for Medical Research, Chapman & Hall, London):

d.intake <- c(5260, 5470, 5640,
  6180, 6390, 6515,
  6805, 7515, 7515,
  8230, 8770)

Необходимо выяснить, отличается ли потребление энергии в группе обследованных женщин от рекомендованного значения 7725 кДж/сутки. Для выполнения теста Уилкоксона в системе R используется функция wilcox.test():

wilcox.test(d.intake, mu = 7725)
 
        Wilcoxon signed rank test with continuity correction
 
data:  d.intake 
V = 8, p-value = 0.0293
alternative hypothesis: true location is not equal to 7725 
 
Warning message:
In wilcox.test.default(d.intake, mu = 7725) :
  cannot compute exact p-value with ties

Видим, что рассчитанное значение критерия V = 8. Вероятность получить такое (или превышающее его) значение при условии, что нулевая гипотеза верна, не превышает 0.05 (p-value = 0.0293). Это позволяет нам отклонить нулевую гипотезу о том, что суточное потребление энергии у обследованных 11 женщин не отличается от принятой нормы. Обратите внимание на выданное программой предупреждение о том, что полученное значени вероятности Р не является точным из-за наличия в данных повторяющихся значений (Warning message... cannot compute exact p-value with ties). Проблема рассчета точных Р-значений при наличии повторяющихся наблюдений в данных характерна для статистических методов, основанных на рангах, и критерий Уилкоксона здесь, увы, не исключение (но см. функцию wilcox_test() из пакета coin). При наличии повторяющихся наблюдений Р-значение рассчитывается путем аппроксимации распредения критерия Уилкоксона нормальным распределением (см. справочный файл по функции ?wilcox.test).


Сравнение двух зависимых выборок

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

Используем рассмотренный ранее пример о суточном потреблении энергии, измеренном у одних и тех же 11 женщин до и после периода менструаций:

data(intake) # из пакета ISwR
attach(intake)
intake
    pre post
1  5260 3910
2  5470 4220
3  5640 3885
4  6180 5160
5  6390 5645
6  6515 4680
7  6805 5265
8  7515 5975
9  7515 6790
10 8230 6900
11 8770 7335

Сравнить два периода по потреблению энергии при помощи критерия Уилкоксона можно следующим образом (обратите внимание на использование аргумента paired = TRUE):

wilcox.test(pre, post, paired = TRUE)
 
 Wilcoxon signed rank test with continuity correction
 
data:  pre and post 
V = 66, p-value = 0.00384
alternative hypothesis: true location shift is not equal to 0 
 
Warning message:
In wilcox.test.default(pre, post, paired = T) :
  cannot compute exact p-value with ties

Как видим, рассчитанное программой P-значение оказалось меньше 0.05, что позволяет нам сделать заключение о наличии статистически значимой разницы в потреблении энергии у исследованных женщин до и после менструации. (Для сравнения: Р-значение, полученное при помощи критерия Стьюдента было << 0.001). Мы можем оценить доверительный интервал, в котором с определенной вероятностью находится истинная величина эффекта, воспользовавшись аргументом conf.int (вероятность задается при помощи аргумента conf.level; по умолчанию рассчитывается 95%-ный доверительный интервал):

wilcox.test(pre, post, paired = TRUE, conf.int = TRUE)
 
 Wilcoxon signed rank test with continuity correction
 
data:  pre and post 
V = 66, p-value = 0.00384
alternative hypothesis: true location shift is not equal to 0 
95 percent confidence interval:
 1037.5 1582.5 
sample estimates:
(pseudo)median 
      1341.332 
 
Warning messages:
1: In wilcox.test.default(pre, post, paired = TRUE, conf.int = TRUE) :
  cannot compute exact p-value with ties
2: In wilcox.test.default(pre, post, paired = TRUE, conf.int = TRUE) :
  cannot compute exact confidence interval with ties

Видим, что истинная разница уровней потребленной энергии с вероятностью 95% находится в интервале от 1037.5 до 1581.5 кДж/сутки (для сравнения: при использовании критерия Стьюдента этот интервал составил 1074.1 - 1566.8 кДж/сутки). Опять-таки, из-за наличия повторяющихся наблюдений, рассчет точных доверительных пределов оказался невозможным (см. пункт 2 в списке предупреждений Warning messages). Псевдомедиана ((pseudo)median) индивидуальных разниц между парными значениями потребления энергии была оценена в 1341.3 кДж/сутки. (Псевдомедианой распределения F называют медиану распределения (u + v) / 2, где u и v являются независимыми переменными, каждая из которых имеет распределение F. Если распределение F симметрично, псевдомедиана и медиана совпадают. Подробнее см. ?wilcox.test).


Сравнение двух независимых выборок

Если сравниваемые выборки являются независимыми (paired = FALSE), то мы имеем дело с критерием Уилкоксона, который в англоязычной литературе называют Wilcoxon rank sum test (как было отмечено выше, этот метод называют также методом Манна-Уитни). Проверяемая с его помощью нулевая гипотеза состоит в том, что центры распределений, из которых происходят сравниваемые выборки, смещены относительно друг друга на величину µ (например, µ = 0). Алгоритм метода состоит в следующем:
  • Все имеющиеся значения ранжируют, игнорируя их знак;
  • Ранги значений, принадлежащих к первой группе, суммируют, получая величину W;
  • Полученное W сравнивают со значением, которое можно было бы ожидать при верной нулевой гипотезе и имеющемся числе степеней свободы; альтернативный подход - находят вероятность случайным образом получить значение W, равное или превышающее наблюдаемое значение (при условии истинности нулевой гипотезы).
Используем рассмотренный раннее пример о суточном расходе энергии (expend) у худощавых женщин (lean) и женщих с избыточным весом (obese):
data(energy) # из пакета ISwR
attach(energy)
str(energy)
'data.frame': 22 obs. of  2 variables:
 $ expend : num  9.21 7.53 7.48 8.08 8.09 ...
 $ stature: Factor w/ 2 levels "lean","obese": 2 1 1 1 1 1 1 1 1 1 ...

Проверим гипотезу об отсутствии разницы в потреблении энергии у женщих из этих двух групп при помощи критерия Уилкоксона для независимых выборок:

wilcox.test(expend ~ stature, paired = FALSE)
 
 Wilcoxon rank sum test with continuity correction
 
data:  expend by stature 
W = 12, p-value = 0.002122
alternative hypothesis: true location shift is not equal to 0 
 
Warning message:
In wilcox.test.default(x = c(7.53, 7.48, 8.08, 8.09, 10.15, 8.4,  :
  cannot compute exact p-value with ties

Согласно полученному Р-значению (p-value = 0.002122), потребление энергии у женщин из рассамтриваемых весовых групп статистически значимо различается. Отвергая нулевую гипотезу о равенстве потребляемой энергии, мы рискуем ошибиться с вероятностью лишь около 0.2%.

Важно отметить одно из ограничений критерия Уилкоксона для двух выборок (зависимых или независимых): если общее количество наблюдений не превышает 6, то обнаружить разницу между выборками с уровнем ошибки в 5% просто невозможно (Dalgaard 2008).


16.05.2012

Вышла новая версия RStudio - v0.96

Два дня тому назад в официальном блоге RStudio было сообщено о выходе новой версии программы. В данном обновлении разработчики одной из самых популярных IDE для R сосредоточились на опциях, направленных на 1) обеспечение воспроизводимости научных исследований путем "грамотного программирования" (соответственно, расширились возможности по работе с Sweave и появилась интеграция с пакетом knitr); 2) создание динамических веб-отчетов в форматах R Markdown и R HTML.

Появилось также несколько новых опций по редактированию кода, включая давно ожидаемое сворачивание кода. Более подробно новшества RStudio описаны в представленном ниже видеоролике (англ. яз.).




09.05.2012

Воспроизводимость результатов при использовании генератора случайных чисел

В практике статистического анализа данных часто приходится иметь дело с необходимостью генерации случайных чисел, подчинящихся тому или иному закону распределения вероятностей (например, при необходимости случайным образом отобрать небольшую выборку из массивной таблицы данных, при использовании бутстреп-методов, методов Монте-Карло, и т.п.).

В системе R реализован целый ряд алгоритмов, позволящих генерировать случайные числа (см. ?RNGkind; по умолчанию используется один из наиболее быстрых генераторов - т.н. "Вихрь Мерсенна" (Mersenne twister)). Очень важно понимать, что практически все эти алгоритмы детерминированы, и поэтому генерируемые с их помощью числа правильно называть "псевдослучайными".

Генератор псевдослучайных чисел (ГПСЧ) начинает свою работу с определенной точки в пространстве возможных чисел. Эта точка называется начальным числом (англ. seed). В R имеется возможность зафиксировать это число так, что при повторном использовании ГПСЧ будет генерироваться точно та же последовательность чисел, что и в первый раз. Это может оказаться полезным в случаях, когда исследователь (по тем или иным причам) желает иметь точную воспроизводимость результатов, получаемых с задействованием ГПСЧ.

В качестве примера создадим таблицу example с двумя столбцами. В первом столбце будут храниться коды уровней гипотетического фактора Factor (три уровня, A, B, и C). Для каждого из этих уровней сгенерируем (псевдо-)случайным образом по 300 нормально распределенных значений с разными средними и стандартными отклонениями:

example = data.frame(
    Factor = rep(c("A", "B", "C"), each = 300),
    Variable = c(rnorm(300, 5, 2),
                 rnorm(300, 4, 3),
                 rnorm(300, 2, 1))
    )
 
 
tapply(example$Variable, example$Factor, summary)
$A
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
-0.9135  3.5950  5.0190  4.8520  6.2250 10.4900 
 
$B
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 -4.407   2.194   4.226   4.028   5.962  13.750 
 
$C
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 -0.517   1.365   1.997   1.975   2.582   4.458


Изобразим на одном графике кривые плотности вероятности для каждой из трех групп А, В, и С (для построения графика испозован ggplot2 - графический пакет для R):

library(ggplot2)
p <- ggplot(example, aes(x = Variable))
p <- p + geom_density(aes(fill = Factor), alpha = 1/2)




Если бы мы повторили создание таблицы example описанным выше способом еще несколько раз, мы получили бы очень похожие, но все-таки несколько различающиеся распределения для групп А, В и С. На следующем рисунке приведен пример для четырех таких таблиц, созданных при помощи ГПСЧ:


Для того, чтобы каждый раз при использовании ГПСЧ в R получать идентичные последовательности чисел, используется функция set.seed() (от set - задать, установить, и seed - начальное число). Как следует из названия, эта функция фиксирует число, служащее начальной точкой для запуска алгоритма генерации (псевдо-)случайных чисел. В качестве аргумента функции указывают любое целое число (не важно, какое именно). Так, при повторном выполнении следующего кода мы будем всегда получать одинаковые последовательности значений для групп А, В и С:

set.seed(1020)
example = data.frame(
    Factor = rep(c("A", "B", "C"), each = 300),
    Variable = c(rnorm(300, 5, 2),
                 rnorm(300, 4, 3),
                 rnorm(300, 2, 1))
    )

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