19 августа 2012

Критерий хи-квадрат для таблиц сопряженности размером больше 2x2



В одном из предыдущих сообщений я описал принцип работы критерия \(\chi^2\) и привел пример его применения в отношении простейшего случая - анализа частот, сведенных в таблицу сопряженности размером \(2\times 2\). Безусловно, на практике можно столкнуться с данными, имеющими более сложную структуру. В целом, можно говорить о таблицах сопряженности размером \(R\times C\), где \(R\) - количество строк, а \(C\) - количество столбцов в таблице. Число степеней свободы в таких таблицах составляет \(df = (R - 1)(C - 1)\).

Источник: wikipedia.org
Предположим, мы выполнили исследование трех популяций наземных моллюсков, в ходе которого отмечали цвет раковины. Цвет выражался следующими тремя градациями - "светлый" (light), "темный" (dark) и "очень темный" (very dark). Необходимо выяснить, различаются ли частоты встречаемости каждого из вариантов окраски в исследованных популяциях.

Допустим, что учтенные частоты выражались следующими значениями:
light <- c(12, 40, 45)
dark <- c(87, 34, 75)
very.dark <- c(3, 8, 2)

Соберем все три вектора в одну матрицу:

color.data <- matrix(c(light, dark, very.dark), nrow = 3,
             dimnames = list(c("Pop1", "Pop2", "Pop3"),
             c("Light", "Dark", "Very dark")))
 
color.data
     Light Dark Very dark
Pop1    12   87         3
Pop2    40   34         8
Pop3    45   75         2

Дальше достаточно подать эту матрицу на функцию chisq.test():

chisq.test(color.data)
 
 Pearson's Chi-squared test
 
data:  color.data 
X-squared = 43.4337, df = 4, p-value = 8.411e-09
 
Warning message:
In chisq.test(color.data) : Chi-squared approximation may be incorrect

Поскольку полученное P-значение намного меньше 0.05 (p-value = 8.411e-09), у нас есть все основания отклонить нулевую гипотезу об отcутствии разницы между популяциями моллюска по частотам встречаемости разных вариантов окраски раковины.


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

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

А вот такой вопрос по теме блога:
В пяти озерах закидывали невод соответственно (23,56,31,18,16)раз. Рыбка-белоглазка выловилась в (0,5,2,0,1) случаях.
Как правильно проверить гипотезу по ХиХи:
X <- c(0,5,2,0,1) # Белоглазка
P <- c(23/144,61/144,33/144,18/144,17/144)
P[5]=1 - sum(P[1:4])
sum(P)
chisq.test(X, p=P)
= X-squared = 3.5589, df = 4, p-value = 0.469
или
X <- t(matrix(c(0,5,2,0,1, # Рыбка встретилась
23,56,31,18,16), ncol=2)) # Рыбка не встретилась
chisq.test(X,correct=FALSE)
= X-squared = 3.3858, df = 4, p-value = 0.4955

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

total <- c(23,56,31,18,16) # общее кол-во попыток
positive <- c(0,5,2,0,1) # кол-во положительных попыток
negative <- total - positive # кол-во отрицательных попыток

X <- matrix(c(positive, negative), ncol = 2)
chisq.test(X)
# ... p-value = 0.4498

# поскольку число наблюдений в ячейках для позитивных случаев
# очень мало (вплоть до 0), то корректнее будет применить точный тест Фишера:
fisher.test(X)
# ... p-value = 0.5578

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

Вообще-то мой вопрос был чисто методический: "в каких случаях используется вектор вероятностей р в параметрах функции chisq.test() и нельзя ли его использовать при обсчете таблиц 2хС, задавая только одну строку и вектор Р".
А зачем мы рассчитываем Хи-квадрат, можно почитать в разделе 3.2 нашей недописанной еще книжки "Рандомизация, бутстреп и методы Монте-Карло", которая лежит на
http://www.ievbras.ru/ecostat/Kiril/Article/A32/Stare.htm
Кстати, почти все примеры там вычисляются в R и коды скриптов представлены в приложении. Так что книга вполне может позиционироваться как практикум по статистике в R для начинающих (каковыми мы и являемся).

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

Вообще-то вопрос был: "Как правильно проверить гипотезу по ХиХи?". На него и был дан ответ (с учетом приведенного примера, в котором, как я понимаю, стояла задача сравнить озера по частоте вылавливания "рыбки").

Что касается второго, уточненного, вопроса, то тут стоит еще раз обратиться к справочному файлу - ?chisq.test, который нам говорит:

"...If x is a matrix with one row or column, or if x is a vector and y is not given, then a goodness-of-fit test is performed (x is treated as a one-dimensional contingency table). The entries of x must be non-negative integers. In this case, the hypothesis tested is whether the population probabilities equal those in p, or are all equal if p is not given."

Другими словами, если задать только один вектор из наблюдаемых частот, то функция будет "думать", что пользователь хочет выполнить goodness-of-fit test (например, в случае проверки соответствия каких-то частот ожидаемому теоретическому вероятностному распределению). При этом если указать вектор с ожидаемыми генеральными вероятностями p, то выборочные частоты будут сравниваться с этими значениями р. Если вектор с р не указан, то предполагается, что все ожидаемые вероятности равны.

Спасибо за ссылку на Вашу книгу (и за труд над ней)! Добавил в раздел "Библиотека".

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

Мне до сих пор непонятно, с какого перепугу гипотезу об однородности вероятностей в таблицах сопряженности считают эквивалентной гипотезе о независимости переменных А и В (см. Афифи и проч.). В нашем примере основная цель выяснить - влияет ли география на встречаемость особей, а не проверка, равны ли все вероятности в совокупности. Думается, надо проверять гипотезу о равенстве вероятностей именно в парах р1=Р1 ; р2=Р2 ; ... , р5 = Р5, где Рi - ожидаемая вероятность вылова в том самом озере (т.е.вектор Р). Проверять же равенство вероятностей р11 = р12 = р21 = р22 = ... = р55 , по моему, - стремление наводить ненужную дополнительную тень на плетень. Так что , при зрелом размышлении, мне кажется логичнее считать как раз по первому варианту. Кстати, прошу прощения, что в условии ошибочно задал не тот вектор частот закидывания невода, нужно (23,61,33,18,17). Впрочем, все это не слишком относится к языку R.

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

Поясните, пожалуйста, следующее - если вывод об ожидаемых значениях частот делается по совокупности попыток - как будет выглядеть корректный код на R?
пусть пример гипотетических данных похож на предыдущий -
color.data

Light Dark Very dark
Pop1 12 87 3
Pop1 17 66 9
Pop1 9 52 4
Pop2 40 34 8
Pop2 45 75 2

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

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

При наличии нескольких повторных измерений, выполненных на одних и тех же экспериментальных единицах (в приведенном Вами примере - это популяции), критерий хи-квадрат для сравнения частот неприменим. Это объясняется тем, что критерий хи-квадрат предполагает независимость наблюдений, тогда как наблюдения выполненные на одних и тех же экспериментальных единицах с большой вероятностью коррелируют, т.е. не являются независимыми (см. обсуждение, посвященное, например, критерию Стьюдента для сравнения двух зависимых выборок: http://r-analytics.blogspot.de/2012/03/t.html).

Если бы Вы имели дело с таблицей сопряженности 2х2, где бинарный признак учтен два раза в одной и той же популяции, для сравнения подошел бы критерий Мак-Немара (https://en.wikipedia.org/wiki/McNemar%27s_test). В R имеется соответствующая функция: http://stat.ethz.ch/R-manual/R-patched/library/stats/html/mcnemar.test.html
Я приведу примеры использования критерия Мак-Немара в одном из будущих сообщений - подписывайтесь на рассылку, чтобы не пропустить :)

Однако анализ, который Вы хотите выполнить, насколько я понимаю, более сложный - есть несколько популяций, обследованных несколько раз. В этом случае Вы автоматически попадаете во Вселенную т.н. "моделей со смешанными эффектами". Более того, поскольку переменная-отклик бинарна, подходящей будет логистическая модель со смешанными эффектами. Литературы по этого рода моделям - масса (Google в помощь!). Отличная отправная точка здесь: http://glmm.wikidot.com/faq

Пакеты R для расчета моделей со смешанными эффектами: nlme и lme4. Есть и другие, но эти два - классика.

По nlme классическая книга: http://books.google.de/books/about/Mixed_Effects_Models_in_S_and_S_PLUS.html?id=y54QDUTmvDcC&redir_esc=y

По lme4 главы недописанной книги можно скачать здесь: http://lme4.r-forge.r-project.org/book/

Удачи!

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

Благодарю Вас, копаю... от забора и в указанном направлении. Действительно, часть данных укладывается в определение "несколько популяций, обследованных несколько раз". Позволю себе уточняющий вопрос - ради понимания содержания группы вариантов анализа, построенных на хи-квадрат (напр. http://www.hr-portal.ru/statistica/gl15/gl15.php), или шкалирование типа [dist.dudi(*.coa)]. Если в нашем модельном случае Pop1 и Pop2 - не отклик, а маркер воздействия (отклик же - собственно распределение частот встречаемости)- а каждое измерение выполнено(примем) на независимой экспериментальной единице - разве не логично попытаться опереться на этот критерий? Ситуация не особо редкая - рассматриваем популяции моллюсков в горной(Pop1) и равнинной(Pop2) зонах, обследовали (однократно) три речки в горной, две в равнинной... при условии, что гипотеза формулируется так же, как выше.
Прошу прощения, если что-то не так понял - пытаюсь разобраться.

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

Если Вы внимательно посмотрите на таблицу сопряженности в примере, приведенном на указанном Вами сайте (http://www.hr-portal.ru/statistica/gl15/gl15.php), то увидите, что степень проявления признака "курение" изучается для нескольких НЕЗАВИСИМЫХ экспериментальных единиц - в данном случае эти единицы представлены группами "старшие менеджеры", "младшие менеджеры" и т.д. Каждая из этих групп людей исследована на предмет курения только ОДИН раз. Поэтому хи-квадрат критерий подходит для анализа таких данных.

Теперь про зоны/речки. Естественно, Pop1 и Pop2 - это не отклик. Отклик, как Вы првильно, заметили - это частота встречаемости. В этом примере Pop1 и Pop2 - это наши экспериментальные единицы, для КАЖДОЙ из которых мы измеряем частоту встречаемости НЕСКОЛЬКО раз - в первом случае три раза (три реки), во втором - 2 раза (две реки). Отсюда и вытекат "проблема" - измерения, выполненные на одних и тех экспериментальных единицах, будут коррелировать и, как следствие, хи-квадрат неприменим. Один из возможных выходов - объединить все данные в пределах одной популяции (т.е. по трем рекам для Pop1 и по двум рекам для Pop2). Тогда можно будет применить хи-квадрат. Такой подход аналогичен усреднению нескольких измерений количественного признака, выполненных на одних и тех же единицах и дальнейшей работе уже с полученными средними. В этом случае теряется часть информации, но зато снимается проблема коррелирующих измерений, выполненных на одних и тех же объектах.

Надеюсь, приведенное объяснение поможет разобраться с Вашей проблемой.

Совет 1: всегда старайтесь иметь четкое представление о том, на каких единицах выполняется измерение того или иного признака.

Совет 2: почитайте весьма полезную классическую статью (> 5000 цитирований!) Халберта про "pseudoreplicates" (псевдоповторы) - http://www.masterenbiodiversidad.org/docs/asig3/Hurlbert_1984_Pseudoreplication.pdf

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

Благодарю, ушёл думать. С работами Хёлберта/Козлова/Петрова знаком, и на этот предмет (случай наблюдаемых/ожидаемых частот) уже просматривал, но явных указаний не нашёл, неявные - сомнительны. Хотелось пройти по кромке корректности - и информацию (о различиях вариаций) не потерять, и с единицами не напутать - да уровень пока не позволяет. Ещё раз спасибо.

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