Методы множественных проверок гипотез, реализованные в пакете multcomp

Последние несколько сообщений были посвящены проблеме множественных проверок статистических гипотез. Для решения этой проблемы разработано большое число методов, различающихся по мощности и применимости в разных ситуациях (так, были рассмотрены методы Бонферрони и Холма, Тьюки, Беньямини-Хохберга и Беньямини-Йекутили). Разнообразие этих методов может создать ощущение неразберихи и привести в замешательство даже опытных исследователей. Тем не менее, между многими методами существует большое сходство. Более того, можно показать, что некоторые методы, известные и используемые под разными названиями и для разных целей, с математической точки зрения эквиваленты (например, тесты Тьюки и Даннета). Используя теорию общих линейных моделей, проф. Франк Брeтц и соавт. (Bretz et al. 2010) разработали общую методологическую схему, объединяющую большинство классических критериев для множественной проверки гипотез. Как это часто происходит в наши дни, соответствующие методологические подходы были реализованы в дополнительном пакете для R - multcomp (от "multiple comparisons" - "множественные сравнения"). Цель данного сообщения - дать описание основных возможностей этого пакета.  Следует подчеркнуть, что это описание будет иметь лишь поверхностный характер. Для полноты картины следует раздобыть указанную выше книгу (Bretz et al. 2010)  - интересующиеся читатели найдут в ней подробные математические выкладки и множество примеров R кода.

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

В качестве примера воспользуемся рассмотренными ранее данными по содержанию стронция в воде пяти водоемов США:

waterbodies <- data.frame(Water = rep(c("Grayson", "Beaver",
                                       "Angler", "Appletree",
                                       "Rock"), each = 6),
                          Sr = c(28.2, 33.2, 36.4, 34.6, 29.1, 31.0,
                                 39.6, 40.8, 37.9, 37.1, 43.6, 42.4,
                                 46.3, 42.1, 43.5, 48.8, 43.7, 40.1,
                                 41.0, 44.1, 46.4, 40.2, 38.6, 36.3,
                                 56.3, 54.1, 59.4, 62.7, 60.0, 57.3)
                          )

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

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

\[y_{i} = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \beta_3 x_3 + \beta_4 x_4 + \epsilon_{i},\]

где \(y_i\) - это i-е значение содержания стронция в воде; \(\beta_0\) - среднее значение содержания стронция для базового уровня изучаемого фактора (в данном случае выбор базового уровня может быть совершенно произвольным; по умолчанию R выберет в качестве базового тот уровень, название которого идет первым по алфавиту - Angler); \(\beta_1\) ... \(\beta_4\) - коэффициенты, отражающие разницу между средним значением содержания стронция в воде "базового водоема" и средними значениями в других водоемах; \(\epsilon_i\) - нормально распределенные остатки.

Параметры приведенной выше модели легко оценить при помощи функции lm() (подробнее см. здесь):

M  <- lm(Sr ~ Water, data = waterbodies)
 
summary(M)
 
Call:
lm(formula = Sr ~ Water, data = waterbodies)
 
Residuals:
    Min      1Q  Median      3Q     Max 
-4.8000 -2.2500 -0.4833  2.2042  5.3000 
 
Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
(Intercept)      44.083      1.276  34.555  < 2e-16 ***
WaterAppletree   -2.983      1.804  -1.654   0.1107    
WaterBeaver      -3.850      1.804  -2.134   0.0428 *  
WaterGrayson    -12.000      1.804  -6.651 5.72e-07 ***
WaterRock        14.217      1.804   7.880 3.09e-08 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
 
Residual standard error: 3.125 on 25 degrees of freedom
Multiple R-squared:  0.8998, Adjusted R-squared:  0.8838 
F-statistic: 56.15 on 4 and 25 DF,  p-value: 3.948e-12

Как видим, в целом средние уровни содержания стронция в воде исследованных водоемов существенно значимо различаются (см. последнюю строку результатов: Р-значение для F-критерия оказалось гораздо меньше 0.05, а именно 3.948е-12). Кроме того, из полученных результатов мы можем сделать определенные выводы касательно того, какие из исследованных водоемов отличались от водоема Angler, автоматически выбранного программой в качестве "эталонного" для сравнений. Так, есть основания считать, что параметры \(\beta_2\), \(\beta_3\) и \(\beta_4\) статистически значимо (на уровне 0.05) отличаются от 0. Другими словами, мы можем утверждать, что среднее содержание стронция в "эталонном" водоеме Angler значимо отличалось от таковых в Beaver (Р = 0.0428), Grayson (P = 5.72e-07) и Rock (P = 3.09e-08).

К сожалению, описанная интерпретация Р-значений, полученных для параметров модели M, имеет один недостаток: эти Р-значения являются безусловными (англ. marginal) (сравните с "условным распределением вероятностей"; англ. conditional probability distribution). Такое название связано с тем, что эти Р-значения рассчитываются исходя из допущения об отсутствии какой-либо связи между параметрами модели. Иными словами, параметры модели рассматриваются как некоррелирущие случайные величины, каждый имеющий свое собственное распределение вероятностей. Проблема, однако, заключается в том, что при оценке параметров модели мы используем один конкретный набор данных. Соответственно, мы не можем строго утверждать, что получаемые параметры не коррелируют друг с другом. Как следствие, проверяя одновременно несколько статистических гипотез (в отношении коррелирующих параметров модели) на одних и тех же данных, мы не выполняем должного контроля над групповой вероятностью ошибки первого рода. Это, в свою очередь, может привести к повышенной вероятности того, что выводы, основанные на таком анализе, не подтвердятся в будущем при сборе дополнительных данных (например, при повторении эксперимента). В ряде областей слабая воспроизводимость результатов может иметь серьезные практические последствия (например, в фармацевтической промышленности, когда какая-либо компания объявляет о повышенной эффективности своего нового дорогого препарата по сравнению с эффективностью более дешевого старого аналога).

Одним из возможных решений указанной проблемы является корректировка полученных Р-значений при помощи того или иного метода множественных сравнений (например, рассмотренная ранее поправка Бонферрони). Однако, как отмечают Bretz et al. (2010), более мощным подходом (в смысле статистической мощности) был бы расчет Р-значений, основанный на непосредственном учете корреляции между параметрами модели. Этот подход как раз и реализован в пакете multcomp. В частности, в зависимости от класса рассматриваемой модели, принимается допущение о том, что наблюдаемые значения ее параметров являются случайными реализациями значений из определенного многомерного распределения (t-распределения или нормального распределения), ковариационная матрица которого отражает корреляцию между параметрами модели. Размерность распределения равна числу одновременно проверяемых статистических гипотез в отношении параметров модели. Соответствующие Р-значения рассчитываются, исходя из свойств этого распределения.

Функция glht()

Функция glht() (от general linear hypothesis testing) из пакета multcomp предоставляет удобный интерфейс для выполнения одновременной проверки статистических гипотез в отношении параметров самых разнообразных статистических моделей, включая общие линейные модели, обобщенные линейные модели, модели со смешанными эффектами и т.д. Основное требование: соответствующий R-объект с результатами расчета той или иной модели должен содержать оценки ее параметров и ковариационную матрицу, которые могут быть извлечены при помощи методов coef и vcov соответственно. Так, для рассчитанной нами выше модели М, имеем:

coef(M)
(Intercept) WaterAppletree    WaterBeaver   WaterGrayson      WaterRock 
      44.08          -2.98          -3.85         -12.00          14.22
 
vcov(M)
               (Intercept) WaterAppletree WaterBeaver WaterGrayson WaterRock
(Intercept)           1.63          -1.63       -1.63        -1.63     -1.63
WaterAppletree       -1.63           3.26        1.63         1.63      1.63
WaterBeaver          -1.63           1.63        3.26         1.63      1.63
WaterGrayson         -1.63           1.63        1.63         3.26      1.63
WaterRock            -1.63           1.63        1.63         1.63      3.26

Синтаксис команд с использованием функции glht() в большинстве случаев имеет следующую структуру:

glht(model,
     linfct,
     alternative = c("two.sided", "less", "greater"),
     ...
)

Здесь аргумент model - это линейная модель, подогнанная с использованием, например, таких стандартных функций, как aov(), lm(), или glm(). На аргумент linfct (от linear function, "линейная функция") подается матрица контрастов, при помощи которой задаются подлежащие проверке гипотезы. Понятие линейных контрастов было подробно рассмотрено ранее. Имеется несколько способов спецификации аргумента linfct - они будут рассмотрены ниже. Аргумент alternative служит для указания типа альтернативных гипотез. Возможные его значения: "two.sided" (двусторонняя), "less" ("меньше чем") и "greater" ("больше чем") (подробнее см. здесь). Многоточие "..." означает, что функция glht() может принимать и другие дополнительные аргументы (см. справочный файл, доступный по команде ?glht).


Примеры использования функции glht()

Продолжая пример с содержанием стронция в водоемах, рассмотрим способы спецификации линейных контрастов, подаваемых на аргумент linfct. Напомню еще раз, что в случае с дисперсионным анализом матрица контрастов используется для кодирования сочетаний экспериментальных групп, подлежащих сравнению (пример матрицы контрастов, заданной пользователем, был приведен ранее). Следует подчеркнуть, что с подлежащими сравнению группами следует  определиться до выполнения анализа данных. Иными словами, проверяемая научная гипотеза (или гипотезы) должна быть сформулирована до сбора данных. Любые гипотезы, формулируемые во время анализа уже имеющихся данных (в англоязычной литературе это называют "data dredging" или "data fishing"), будут диктоваться свойствами именно этих данных, и нет никакой гарантии, что при повторении эксперимента исследователь обнаружит те же закономерности (это распространенная проблема, известная как "overfitting", или "переобучение модели").

В случае с дисперсионным анализом, самый простой способ задать матрицу контрастов состоит в использовании служебной функции mcp() (от "multiple comparisons"), которая сделает всю работу автоматически. От нас требуется лишь указать фактор, уровни которого соответствуют сравниваемым группам, и название одного из "встроенных" критериев для выполнения сравнений. Так, для выполнения всех попарных сравнений средних значений концентрации стронция в исследованных водоемах при помощи критерия Тьюки, соответствующая команда будет иметь вид:

glht(M, linfct = mcp(Water = "Tukey"))
 
  General Linear Hypotheses
 
Multiple Comparisons of Means: Tukey Contrasts
 
 
Linear Hypotheses:
                         Estimate
Appletree - Angler == 0   -2.9833
Beaver - Angler == 0      -3.8500
Grayson - Angler == 0    -12.0000
Rock - Angler == 0        14.2167
Beaver - Appletree == 0   -0.8667
Grayson - Appletree == 0  -9.0167
Rock - Appletree == 0     17.2000
Grayson - Beaver == 0     -8.1500
Rock - Beaver == 0        18.0667
Rock - Grayson == 0       26.2167

Результатом выполнения приведенной команды является таблица с перечисленными проверяемыми гипотезами (например, первая гипотеза состоит в том, что средние значения концентрации стронция в водоемах Appletree и Angler не различаются: Appletree - Angler == 0) и наблюдаемыми уровнями различий между группами (так, в водоеме Appletree содержание стронция в среднем оказалось на 2.983 мг/мл меньше, чем в Angler). Позднее мы увидим, как можно получить также Р-значения для каждой из этих гипотез. Помимо критерия Тьюки, реализованы и другие широко используемые критерии, как например, критерий Даннетта ("Dunnett") (применяется для сравнения всех групп с контрольной) и критерий Уильямса ("Williams") (подобно критерию Даннетта, применяется для сравнения всех групп с контрольной, но только в случае, если групповые средние демонстрируют монотонный тренд на убывание или возрастание; широко используется в токсикологии при анализе зависимости величины биологического эффекта от дозы исследуемого вещества). С полным списком "встроенных" критериев можно ознакомиться в справочном файле, доступном по команде ?contrMat.

При необходимости выполнить сравнения только для определенных групп, мы можем воспользоваться другим способом спецификации матрицы контрастов. Этот способ состоит в подаче на функцию mcp() списка из проверяемых гипотез, используя символьные выражения (англ. expression). Предположим, нас интересуют сравнения только следующих водоемов друг с другом: "Rock vs. Angler", "Grayson vs. Appletree" и "Grayson vs. Beaver". Соответствующая команда будет иметь вид:

glht(M, linfct = mcp(Water = c(
  "Rock - Angler = 0",
  "Grayson - Appletree = 0",
  "Grayson - Beaver = 0")
  )
)
 
  General Linear Hypotheses
 
Multiple Comparisons of Means: User-defined Contrasts
 
 
Linear Hypotheses:
                         Estimate
Rock - Angler == 0         14.217
Grayson - Appletree == 0   -9.017
Grayson - Beaver == 0      -8.150

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

contr <- rbind("Rock - Angler" = c(-1, 0, 0, 0, 1),
               "Grayson - Appletree" = c(0, -1, 0, 1, 0),
               "Grayson - Beaver" = c(0, 0, -1, 1, 0)
               )
 
contr
                   [,1] [,2] [,3] [,4] [,5]
Rock - Angler         -1    0    0    0    1
Grayson - Appletree    0   -1    0    1    0
Grayson - Beaver       0    0   -1    1    0

Теперь подадим матрицу contr на функцию mcp():

glht(M, linfct = mcp(Water = contr))
 
  General Linear Hypotheses
 
Multiple Comparisons of Means: User-defined Contrasts
 
 
Linear Hypotheses:
                         Estimate
Rock - Angler == 0         14.217
Grayson - Appletree == 0   -9.017
Grayson - Beaver == 0      -8.150

Как видим, результат оказался идентичным предыдущему.


Метод summary

Метод summary позволяет получить дополнительные результаты расчетов, выполняемых функцией glht(). В частности, мы можем вывести Р-значения для соответствующих сравнений. Например:

summary(glht(M, linfct = mcp(Water = "Tukey")))
 
  Simultaneous Tests for General Linear Hypotheses
 
Multiple Comparisons of Means: Tukey Contrasts
 
 
Fit: lm(formula = Sr ~ Water, data = waterbodies)
 
Linear Hypotheses:
                         Estimate Std. Error t value Pr(>|t|)    
Appletree - Angler == 0   -2.9833     1.8042  -1.654  0.47913    
Beaver - Angler == 0      -3.8500     1.8042  -2.134  0.23757    
Grayson - Angler == 0    -12.0000     1.8042  -6.651  < 1e-04 ***
Rock - Angler == 0        14.2167     1.8042   7.880  < 1e-04 ***
Beaver - Appletree == 0   -0.8667     1.8042  -0.480  0.98848    
Grayson - Appletree == 0  -9.0167     1.8042  -4.998  0.00035 ***
Rock - Appletree == 0     17.2000     1.8042   9.533  < 1e-04 ***
Grayson - Beaver == 0     -8.1500     1.8042  -4.517  0.00112 ** 
Rock - Beaver == 0        18.0667     1.8042  10.014  < 1e-04 ***
Rock - Grayson == 0       26.2167     1.8042  14.531  < 1e-04 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- single-step method)

Кроме получения подобной сводной информации, мы можем выполнить и некоторые другие виды анализа, воспользовавшись аргументом test функции summary (подробнее см. справочный файл - ?summary.glht).


Метод confint

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

mult <- glht(M, linfct = mcp(Water = contr))
 
confint(mult, level = 0.95)
 
  Simultaneous Confidence Intervals
 
Multiple Comparisons of Means: User-defined Contrasts
 
 
Fit: lm(formula = Sr ~ Water, data = waterbodies)
 
Quantile = 2.5324
95% family-wise confidence level
 
 
Linear Hypotheses:
                         Estimate lwr      upr     
Rock - Angler == 0        14.2167   9.6478  18.7855
Grayson - Appletree == 0  -9.0167 -13.5855  -4.4478
Grayson - Beaver == 0     -8.1500 -12.7188  -3.5812

В выведенных результатах столбцы lwr и upr содержат нижний и верхний доверительный пределы соответствующих оценок параметров модели (Estimate).

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

plot(confint(mult, level = 0.95))

На этом я завершаю обзор возможностей пакета multcomp, равно как и тему множественных проверок статистических гипотез. Как было отмечено выше, это сообщение дает лишь краткое введение в возможности multcomp. Детальное изложение представлено в книге Bretz et al. (2010). Кроме того, на странице этого пакета на сайте CRAN можно найти несколько руководств, содержащих большое количество примеров (см. раздел Vignettes).