В одном из предыдущих сообщений я описал, как в R можно рассчитать критерий хи-квадрат для таблиц сопряженности размером 2 х 2. Как правило, таблицу сопряженности 2 х 2 получают в ходе единичного эксперимента (или обсервационного исследования), направленного на изучение распределения того или иного бинарного признака в двух группах объектов (например, в экспериментальной и контрольной группах). Но что, если один и тот же эксперимент повторяют несколько раз? Например, в ходе клинических испытаний часто эффективность какого-либо нового препарата исследуют по одинаковой схеме в разных медицинских учреждениях. В результате получают набор из К таблиц сопряженности размером 2 х 2, где К - это количество участвовавших в исследовании медицинских центров. По разным причинам можно ожидать, что результаты эксперимента будут несколько варьировать от центра к центру. Соответственно, "медицинский центр" становится важной ковариатой, действие которой мы должны учесть при установлении эффективности испытываемого нового препарата. Одним из статистических методов, позволяющих это сделать, является рассмотренный ниже критерий Кохрана-Мантеля-Хензеля (англ. "Cochran-Mantel-Haenszel test" или просто "CMH test" - по фамилии авторов Cochran (1954) и Mantel and Haenszel (1959)). (Насколько мне известно, устоявшегося перевода названия этого критерия в русскоязычной литературе нет - кроме приведенного, встречаются, например, названия "критерий Кохрана-Мантеля-Гензеля" и "критерий Кохрана-Мантеля-Хенселя". Обсуждение того, как правильно перевести фалимию последнего автора, можно найти здесь).
Принцип работы СМН-критерия и R-функцию для его вычисления я опишу на примере из книги Алана Агрести (Agresti 2002, стр. 230). Имеются результаты исследования эффективности действия мази, содержащей новое действующее вещество, в отношении определенного инфекционного заболевания. Дизайн экперимента заключался в следующем: две группы испытуемых подвергались лечению при помощи нового препарата и плацебо. В конце эксперимента учтено количество случаев успешного и неуспешного лечения в обоих группах. Эксперимент был повторен по описанной схеме в 8 медицинских центрах, при разном общем количестве испытуемых и разном их количестве в отдельных экспериментальных группах в каждом центре. Данные представлены в таблице ниже:
Медицинский центр
|
Экспериментальная группа
|
Исход лечения
|
|
Успешное
|
Неуспешное
|
||
1
|
Препарат
Контроль |
11
10 |
25
27 |
2
|
Препарат
Контроль |
16
22 |
4
10 |
3
|
Препарат
Контроль |
14
7 |
5
12 |
4
|
Препарат
Контроль |
2
1 |
14
16 |
5
|
Препарат
Контроль |
6
0 |
11
12 |
6
|
Препарат
Контроль |
1
0 |
10
10 |
7
|
Препарат
Контроль |
1
1 |
4
8 |
8
|
Препарат
Контроль |
4
6 |
2
1 |
Обозначим ячейки каждой k-й таблицы сопряженности следующим образом:
Экспериментальная группа
|
Успешное лечение
|
Неуспешное лечение
|
Всего
|
Препарат
|
\(n_{k11}\)
|
\(n_{k12}\)
|
\(n_{k1}.\)
|
Контроль
|
\(n_{k21}\)
|
\(n_{k22}\)
|
\(n_{k2}.\)
|
Всего
|
\(n_{k}._{1}\)
|
\(n_{k}._{2}\)
|
\(n_{k}..\)
|
Тогда формулу CMH-критерия можно записать так (если вы не можете рассмотреть отдельные элементы формулы из-за малого размера, кликните по ней правой кнопкой мыши и воспользуйтесь опцией Math Settings > Zoom Trigger для настройки zoom-режима):
\[CMH = \frac{(|\sum_{k}(n_{k11}-\frac{n_{k}._{1}n_{k1}.}{n_{k}..})| - \frac{1}{2})^2}{\sum_{k}\frac{n_{k}._{1}n_{k}._{2}n_{k1}.n_{k2}.}{n^2_{k}..(n_{k}..-1)}}\]
Как видно из приведенной формулы, СМН-критерий очень схож с критерием хи-квадрат. В частности, в числителе формулы мы видим сумму квадратов разностей между наблюдаемой частотой в ячейке \(n_{k11}\) и ее ожидаемой частотой при верной нулевой гипотезе. Знаменатель содержит оценку дисперсии этих квадратов разностей. В связи с таким сходством, CMH-критерий часто называют еще критерием хи-квадрат Мантеля-Хензеля. Нулевая гипотеза, проверяемая при помощи CMH-критерия, заключается в том, что между двумя анализируемыми качественными признаками (в нашем случае "экспериментальная группа" и "исход лечения") нет никакой связи. (Технически говоря, тестируется гипотеза о том, что отношение шансов в каждой из K таблиц сопряженности равно 1; про отношение шансов см., например, здесь и здесь).
В базовой версии R для выполнения CMH-теста имеется функция mantelhaen.test(). Данные подаются на нее в виде трехмерного массива, состоящего из 2 х 2 х K элементов. В нашем случае данные можно ввести следующим образом:
drug <- array(c(11, 10, 25, 27, 16, 22, 4, 10, 14, 7, 5, 12, 2, 1, 14, 16, 6, 0, 11, 12, 1, 0, 10, 10, 1, 1, 4, 8, 4, 6, 2, 1), dim = c(2, 2, 8), dimnames = list( Group = c("Drug", "Control"), Response = c("Success", "Failure"), Center = c("1", "2", "3", "4", "5", "6", "7", "8"))) drug , , Center = 1 Response Group Success Failure Drug 11 25 Control 10 27 , , Center = 2 Response Group Success Failure Drug 16 4 Control 22 10 , , Center = 3 Response Group Success Failure Drug 14 5 Control 7 12 , , Center = 4 Response Group Success Failure Drug 2 14 Control 1 16 , , Center = 5 Response Group Success Failure Drug 6 11 Control 0 12 , , Center = 6 Response Group Success Failure Drug 1 10 Control 0 10 , , Center = 7 Response Group Success Failure Drug 1 4 Control 1 8 , , Center = 8 Response Group Success Failure Drug 4 2 Control 6 1
Далее подаем массив drug на функцию mantelhaen.test():
mantelhaen.test(drug) Mantel-Haenszel chi-squared test with continuity correction data: drug Mantel-Haenszel X-squared = 5.6716, df = 1, p-value = 0.01724 alternative hypothesis: true common odds ratio is not equal to 1 95 percent confidence interval: 1.177590 3.869174 sample estimates: common odds ratio 2.134549
Как видим, исход лечения оказался статистически значимо связанным с тем, получали ли испытуемые новый препарат (p-value = 0.01724). Иными словами, мы можем принять альтернативную гипотезу, которая заключается в том, что истинное общее отношение шансов (для всех таблиц) не равно 1 (alternative hypothesis: true common odds ratio is not equal to 1) . Это подтверждает также рассчитанный программой 95%-ный доверительный интервал для отношения шансов, который не включает 1 (95 percent confidence interval: 1.177590 - 3.869174).
Графически мы можем изобразить имеющиеся данные следующим образом:
library(reshape) # для функции melt() drug.df <- data.frame( melt(drug, id=c("Center", "Group", "Response")) ) library(ggplot2) # графический пакет p <- ggplot(data = drug.df, aes(x = Center, y = value, fill = Response)) + ylab("%") p + geom_bar(stat = "identity", position = "fill") + facet_grid(Group~.)
При интерпретации результатов CMH-теста, выполняемого при помощи функции mantelhaen.test(), следует иметь в виду следующее:
- Выполняется ли поправка нa непрерывность? (См. формулу критерия выше - поправка выполняется за счет отнятия 0.5 в числителе.) По умолчанию поправка используется. При необходимости ее отключить воспользуйтесь аргументом correction = FALSE (получаемое в результате этого P-значение может измениться, что в некоторых случаях может привести к противоположному заключению относительно справедливости нулевой гипотезы).
- Какой вариант гипотезы проверяется? Возможные варианты: altenative = "two.sided" (двухсторонний критерий; по умолчанию), alternative = "greater" (больше) или alternative = "less" (меньше).
- Однородны ли отношения шансов в каждой из анализируемых таблиц сопряженности? Функции mantelhaen.test() автоматически "предполагает", что это условие выполняется. По поводу необходимости проверки однородности отношений шансов однозначного мнения нет (подробнее см. здесь). Тем не менее, функции cmh_test() из пакета coin и epi.2by2() из пакета epiR позволяют выполнить специально предназначенный для этого тест Бреслоу-Дэя (Breslow and Day 1980).
- Наблюдается ли сходный отклик на экспериментальное воздействие во всех имеющихся таблицах сопряженности? Мощность CMH-критерия значительно снижается, когда наблюдаются разнонапрвленные отклики в анализируемых таблицах сопряженности (например, если бы в нашем случае в части таблиц наблюдалось выраженное снижение, а в части таблиц - выраженное повышение частоты случаев успешного лечения среди испытуемых, получавших новый препарат). Это, соответственно, затрудняет интерпретацию высоких Р-значений (скажем, >0.05): вызвано ли наблюдаемое высокое Р-значение отсутствием эффекта или это результат наличия разнонаправленных откликов?
Литература:
- Agresti A (2002) Categorical data analysis. 2nd ed. John Wiley and Sons, 737 pp.
- Breslow NE, Day NE (1980) Statistical methods in cancer research I: The Analysis of Case-Control Studies. International Agency for Research on Cancer, Lyon
- Cochran WG (1954) Some methods for strengthening the common χ2 tests. Biometrics 10: 417-451
- Mantel N, Haenszel W (1959) Statistical aspects of the analysis of data from retrospective studies of disease. J. Natl. Cancer Inst. 22: 719-748
Отправить комментарий