Как было показано ранее, однофакторный дисперсионный анализ (ANOVA) представляет собой частный случай общей линейной модели, в которой единственный предиктор представлен категориальной переменной (фактором) с несколькими уровнями (2 и более). В случае многофакторного дисперсионного анализа имеется два или более интересующих нас фактора. Категориальные предикторы могут быть также включены в модели с количественными предикторами, и тогда мы будем иметь дело с ковариационным анализом. Важным понятием при работе с категориальными предикторами, которому, к сожалению, уделяется недостаточно внимания в соответствующей методической литературе, является понятие "контрастов" (англ. contrasts). Ниже я постараюсь дать небольшое введение на эту тему и привести примеры применения контрастов в R. Для простоты изложения речь будет идти только об однофакторном дисперсионном анализе.
Определение понятия
Контраст представляет собой взвешенную сумму средних значений каждого уровня какой-либо категориальной переменной, где сумма весовых коэффициентов обычно (но не всегда) равна 0. Так, в выражении
\[L = c_1 \bar{x}_1 + c_2 \bar{x}_2 \dots + c_k \bar{x}_k \]
\(c_1 \dots c_k\) - это весовые коэффициенты для каждого из k уровней изучаемого фактора, а \(\bar{x}_1 \dots \bar{x}_k\) - средние значения по каждому из этих уровней.
В зависимости от того, как заданы их коэффициенты, контрасты позволяют проверять разные гипотезы о средних значениях отдельных уровней или комбинаций уровней анализируемого фактора (или факторов). Как отмечено выше, обычно весовые коэффициенты контраста должны в сумме давать 0. Для этого должны выполняться следующие обязательные условия:
- Коэффициенты сравниваемых уровней (или комбинаций уровней) должны иметь разные знаки (+ или -)
- Коэффициенты уровней, не представляющих для нас интереса, приравниваются нулю.
[,1] A -1 B 0 C 1
Аналогично, при сравнении средних значений уровней А и В коэффициенты соответствующего контраста составят:
[,1] A -1 B 1 C 0
Обратите внимание на то, что выбор значений весовых коэффициентов в контрастах совершенно произволен - главное, чтобы в сумме они давали 0. Так, в первом из приведенных примеров коэффициенты вполне могли бы составить
[,1] A -0.5 B 0.0 C 0.5
Также обратите внимание на то, что при наличии трех уровней фактора, мы можем иметь только два контраста. Это и неудивительно, поскольку сравнив А с В и А с С, мы, пусть и косвенным образом, но уже сравнили и В с С. Поэтому при попарном сравнении средних значений фактора с тремя уровнями мы получаем следующую матрицу контрастов:
При необходимости мы можем сравнивать также комбинации из нескольких уровней фактора с другими его уровнями (или комбинациями уровней), что, как несложно догадаться, открывает широкие возможности для проверки самых разнообразных гипотез в отношении анализируемых данных. Так, по тем или иным причинам, мы могли бы захотеть сравнить среднее значение уровня А с общим средним значением уровней В и С. Тогда коэффициенты контраста можно было бы выразить следующим образом:
Стандартные контрасты в R
В базовой версии R реализованы четыре основных типа контрастов:
[,1] A 2 B -1 C -1
Стандартные контрасты в R
В базовой версии R реализованы четыре основных типа контрастов:
- treatment contrasts: контрасты комбинаций условий
- sum contrasts: контрасты сумм
- Helmert contrasts: контрасты Хелмерта
- polynomial contrasts: полиномиальные контрасты для упорядоченных факторов (англ. ordered factors) (т.е. таких факторов, чьи уровни могут быть естественным образом упорядочены по возрастанию или убыванию; например, в отношении ежемесячного дохода можно выделить уровни "низкий" -> "средний" -> "высокий"). Полиномиальные контрасты применяются гораздо реже других и здесь не рассматриваются. Подробнее об этом, а также о других типах контрастов, можно почитать здесь (англ. яз.).
По умолчанию в R для неупорядоченных факторов применяются "контрасты комбинаций условий", а для упорядоченных факторов - полиномиальные контрасты. Мы можем выяснить текущие глобальные настройки контрастов при помощи функции options():
options()$contrasts unordered ordered "contr.treatment" "contr.poly"
Та же функция позволяет изменить глобальные настройки контрастов:
Контрасты комбинаций условий (treatment contrasts)
Как отмечено выше, этот тип контрастов используется в R по умолчанию. Собственно, о лежащем в его основе принципе я уже писал при рассмотрении дисперсионного анализа как частного случая общей линейной модели. Применение этого типа контрастов предполагает, что один из уровней изучаемого фактора (например, контроль при анализе экспериментальных данных) принимается в качестве базового, а все остальные уровни сравниваются с этим базовым уровнем.
В качестве примера используем данные эксперимента по изучению эффективности 6 разных инсектицидов (таблица InsectSprays с этими данными входит в состав базовой версии R; см. также здесь). Каждым из этих средств обработали по 12 растений, после чего подсчитали количество выживших на растениях насекомых.
InsectSprays count spray 1 10 A 2 7 A 3 20 A 4 14 A 5 14 A 6 12 A 7 10 A 8 23 A 9 17 A 10 20 A 11 14 A 12 13 A 13 11 B 14 17 B 15 21 B ... ... ... 69 26 F 70 26 F 71 24 F 72 13 F
boxplot(count ~ spray, data = InsectSprays, xlab = "Инсектицид", ylab = "Число выживших насекомых", col = "coral")
Мы можем применить дисперсионный анализ для проверки нулевой гипотезы об отсутствии разницы между инсектицидами по их эффективности:
M1 <- aov(count ~ spray, data = InsectSprays) summary(M1) Df Sum Sq Mean Sq F value Pr(>F) spray 5 2669 533.8 34.7 <2e-16 *** Residuals 66 1015 15.4 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Согласно полученному результату (Pr(>|t|) <2e-16), мы делаем заключение о том, исследованные инсектициды существенно различаются по эффективности действия (подробнее об интерпретации выдаваемой программой информации см. здесь). Однако выполненный таким образом дисперсионный анализ не позволяет выяснить, какие именно пары инсектицидов различались между собой. Для этого следует выполнить анализ путем подгонки линейной модели при помощи функции lm() (подробнее см. здесь):
M2 <- lm(count ~ spray, data = InsectSprays) summary(M2) Call: lm(formula = count ~ spray, data = InsectSprays) Residuals: Min 1Q Median 3Q Max -8.333 -1.958 -0.500 1.667 9.333 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 14.5000 1.1322 12.807 < 2e-16 *** sprayB 0.8333 1.6011 0.520 0.604 sprayC -12.4167 1.6011 -7.755 7.27e-11 *** sprayD -9.5833 1.6011 -5.985 9.82e-08 *** sprayE -11.0000 1.6011 -6.870 2.75e-09 *** sprayF 2.1667 1.6011 1.353 0.181 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 3.922 on 66 degrees of freedom Multiple R-squared: 0.7244, Adjusted R-squared: 0.7036 F-statistic: 34.7 on 5 and 66 DF, p-value: < 2.2e-16
Как видим, программа автоматически выбрала в качестве базового уровня группу наблюдений для препарата A (поскольку буква A идет первой по алфавиту) - см. строку, озаглавленную как (Intercept) в таблице с результатами анализа. Среднее число выживших насекомых в этой группе составило 14.5. В следующей строке приведена информация, отражающая, насколько препарат B эффективнее по сравнению с препаратом А: видим, что среднее количество выживших насекомых в группе B было несколько выше, чем в группе A (на 0.83), но это повышение не было статистически значимым (Pr(>|t|) = 0.604). В третьей строке представлена информация о различиях между базовым уровнем (препарат А) и препаратом C: эффективность препарата С существенно (Pr(>|t|) <7.27e-11) выше, чем эффективность препарата А (число выживших насекомых в группе С в среднем оказалось на 12.4 ниже, чем в группе А, т.е. 14.5 - 12.4 = 2.1). Интерпретация следующих строк таблицы с результатами анализа аналогична.
Функция contrasts() позволяет ознакомиться с матрицей контрастов для того или иного фактора:
Это - типичный пример матрицы, содержащей весовые коэффициенты "контрастов комбинаций условий". Как отмечено выше, суммы весовых коэффициентов по каждому столбцу матрицы контрастов не всегда равны нулю. Контрасты комбинаций условий - как раз тот случай, поскольку все сравнения выполняются относительно базового уровня (препарат А в рассматриваемом примере).
Следующая таблица обобщает принцип, заложенный в основу контрастов комбинаций условий (использованные в таблице обозначения обсуждались ранее здесь):
Параметр модели
|
Оценка
|
Проверяемая нулевая гипотеза
|
\(\beta_0\)
|
среднее значение для базового уровня: \(\bar{x}_1\)
|
\(H_0: \mu_1 = 0\)
|
\(\beta_1\)
|
среднее значение второго уровня за вычетом среднего
значения базового уровня: \(\bar{x}_2 - \bar{x}_1\)
|
\(H_0: \mu_2 - \mu_1 = 0\)
|
\(\beta_2\)
|
среднее значение третьего уровня за вычетом среднего
значения базового уровня: \(\bar{x}_3 - \bar{x}_1\)
|
\(H_0: \mu_3 - \mu_1 = 0\)
|
…
|
…
|
…
|
…
|
…
|
…
|
Контрасты сумм (sum contrasts)
При параметризации модели дисперсионного анализа с помощью контрастов сумм базовый уровень, с которым сравниваются остальные уровни, представляет собой среднее значение из средних по каждой группе. Разберемся с этим типом контрастов, продолжив пример с инсектицидами.
Используя функцию constrasts() в сочетании с contr.sum() (подробнее см. ?constr.sum) мы можем изменить исходное поведение R в плане того, какие контрасты программа автоматически применяет в отношении конкретного фактора (обратите внимание: для изменения глобальных настроек следует использовать функцию options() - см. выше):
contrasts(InsectSprays$spray) <- contr.sum(n = 6)
В итоге получаем матрицу контрастов сумм:
contrasts(InsectSprays$spray) [,1] [,2] [,3] [,4] [,5] A 1 0 0 0 0 B 0 1 0 0 0 C 0 0 1 0 0 D 0 0 0 1 0 E 0 0 0 0 1 F -1 -1 -1 -1 -1
Как видим, в отличие от матрицы с контрастами комбинаций условий, суммы по всем столбцам матрицы контрастов сумм равны нулю.
Выполнив дисперсионный анализ, увидим, что интерпретация оцененных параметров модели теперь также изменилась:
M3 <- lm(count ~ spray, data = InsectSprays) summary(M3) Call: lm(formula = count ~ spray, data = InsectSprays) Residuals: Min 1Q Median 3Q Max -8.333 -1.958 -0.500 1.667 9.333 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 9.5000 0.4622 20.554 < 2e-16 *** spray1 5.0000 1.0335 4.838 8.22e-06 *** spray2 5.8333 1.0335 5.644 3.78e-07 *** spray3 -7.4167 1.0335 -7.176 7.87e-10 *** spray4 -4.5833 1.0335 -4.435 3.57e-05 *** spray5 -6.0000 1.0335 -5.805 2.00e-07 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 3.922 on 66 degrees of freedom Multiple R-squared: 0.7244, Adjusted R-squared: 0.7036 F-statistic: 34.7 on 5 and 66 DF, p-value: < 2.2e-16
Первая строка в таблице с результатами анализа - (Intercept) - содержит среднее значение числа выживших насекомых, рассчитанное по средними значениям каждой группы. Мы можем легко это проверить, выполнив команду
with(InsectSprays, mean(tapply(count, spray, mean))) [1] 9.5
Во второй строке содержится информации относительно того, насколько среднее число выживших насекомых в следующей группе (B) отличается от среднего значения, вычисленного для всех групп ("общее среднее"). Аналогично, в третьей строке содержится информации о том, насколько среднее число выживших насекомых в группе С отличается от общего среднего, и т.д. Заметьте: по причинам, описанным выше в разделе "Определение понятия", в таблице с результатами анализа представлено только пять сравнений с общим средним значением.
Следующая таблица обобщает принцип, заложенный в основу контрастов сумм:
Параметр модели
|
Оценка
|
Проверяемая нулевая гипотеза
|
\(\beta_0\)
|
среднее значение из средних по каждой группе ("общее среднее"): \(\sum{\bar{x}_i} / k\)
| \(H_0: \sum{\bar{x}_i} / k = 0\) |
\(\beta_1\)
|
среднее значение первого уровня за вычетом общего среднего: \(\bar{x}_1 - \sum{\bar{x}_i} / k\)
| \(H_0: \bar{x}_1 - \sum{\bar{x}_i} / k = 0\) |
\(\beta_2\)
|
среднее значение второго уровня за вычетом среднего значения общего среднего: \(\bar{x}_2 - \sum{\bar{x}_i} / k\)
| \(H_0: \bar{x}_2 - \sum{\bar{x}_i} / k = 0\) |
…
|
…
|
…
|
…
|
…
|
…
|
Контрасты Хелмерта
При параметризации модели дисперсионного анализа с помощью контрастов Хелмерта базовый уровень представляет собой среднее значение из средних по каждой группе ("общее среднее"). Сравнения групповых средних с этим базовым уровнем выполняются согласно схеме, обобщенной в следующей таблице:
Параметр модели
|
Оценка
|
Проверяемая нулевая гипотеза
|
\(\beta_0\)
|
среднее значение из средних по каждой группе ("общее среднее"): \(\sum{\bar{x}_i} / k\)
| \(H_0: \sum{\bar{x}_i} / k = 0\) |
\(\beta_1\)
| \(\bar{x}_2 - (\bar{x}_1 + \sum{\bar{x}_i} / k)\) | \(H_0: \bar{x}_2 - (\bar{x}_1 + \sum{\bar{x}_i} / k) = 0\) |
\(\beta_2\)
| \(\bar{x}_3 - (\bar{x}_1 + \bar{x}_2 + \sum{\bar{x}_i} / k)\) | \(H_0: \bar{x}_3 - (\bar{x}_1 + \bar{x}_2 + \sum{\bar{x}_i} / k) = 0\) |
…
|
…
|
…
|
…
|
…
|
…
|
Применим контрасты Хелмерта к данным по инсектицидам:
# изменение контрастов для фактора spray на контрасты Хелмерта: contrasts(InsectSprays$spray) <- contr.helmert(n = 6) contrasts(InsectSprays$spray) [,1] [,2] [,3] [,4] [,5] A -1 -1 -1 -1 -1 B 1 -1 -1 -1 -1 C 0 2 -1 -1 -1 D 0 0 3 -1 -1 E 0 0 0 4 -1 F 0 0 0 0 5
Обратите внимание: суммы по всем столбцам матрицы контрастов Хелмерта равны нулю. Кроме того, сумма произведений элементов любых двух столбцов матрицы также равна нулю. Например:
mat <- contrasts(InsectSprays$spray) sum(mat[, 1]*mat[, 2]) [1] 0
Контрасты с такими свойствами называются ортогональными.
Рассчитаем теперь параметры новой линейной модели:
M4 <- lm(count ~ spray, data = InsectSprays) summary(M4) Call: lm(formula = count ~ spray, data = InsectSprays) Residuals: Min 1Q Median 3Q Max -8.333 -1.958 -0.500 1.667 9.333 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 9.5000 0.4622 20.554 < 2e-16 *** spray1 0.4167 0.8006 0.520 0.604 spray2 -4.2778 0.4622 -9.255 1.53e-13 *** spray3 -1.4306 0.3268 -4.377 4.38e-05 *** spray4 -1.1417 0.2532 -4.510 2.73e-05 *** spray5 1.4333 0.2067 6.934 2.12e-09 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 3.922 on 66 degrees of freedom Multiple R-squared: 0.7244, Adjusted R-squared: 0.7036 F-statistic: 34.7 on 5 and 66 DF, p-value: < 2.2e-16
Интрепретация полученных параметров модели выполняется в соответствии с последней из приведенных выше таблиц. Следует отметить, что несмотря на свои привлекательные математические свойства (ортогональность), контрасты Хелмерта редко применяются на практике поскольку получаемые параметры модели сложно интерпретировать в приложении к конкретным проблемам. Тем не менее, в некоторых статистических программах этот тип контрастов является настройкой по умолчанию - например, в S-PLUS - коммерческой реализации языка S, который является предшественником R.
Контрасты, заданные пользователем
Помимо рассмотренных стандартных типов контрастов, пользователь имеет возможность задавать свои собственные контрасты, предназначенные для проверки конкретных гипотез в отношении анализируемых данных.
Продолжая работать с данными по эффективности инсектицидов, предположим, что нас интересуют ответы на три конкретных вопроса (пример заимствован отсюда):
- Различаются ли исследованные препараты по своей эффективности в целом? (Выше мы уже выяснили, что так оно и есть, но допустим, что пока мы этого не знаем.)
- Различается ли средняя эффективность трех первых инсектицидов (А, В и С) от трех последних (D, E и F)?
- Имеются основания ожидать, что средняя эффективность препаратов А, В и F в среднем ниже (т.е. количество выживших насекомых после обработки ими выше), чем эффективность трех остальных препаратов. Так ли это?
# коэффициенты контраста для ответа на второй вопрос: con1 <- c(1, 1, 1, -1, -1, -1) # коэффициенты контраста для ответа на третий вопрос: con2 <- c(1, 1, -1, -1, -1, 1) # объединение обоих векторов в матрицу контрастов: con.matrix <- cbind(con1, con2) con.matrix con1 con2 [1,] 1 1 [2,] 1 1 [3,] 1 -1 [4,] -1 -1 [5,] -1 -1 [6,] -1 1
Сообщим программе, что созданная нами матрица контрастов con.matrix должна использоваться для всех новых моделей, включающих фактор spray:
contrasts(InsectSprays$spray) <- con.matrix
Выполним дисперсионный анализ при помощи функции aov():
M5 <- aov(count ~ spray, data = InsectSprays)
Как обычно, результаты анализа можно вывести на экран при помощи функции summary(). Однако, учитывая специфичность выполняемых межгрупповых сравнений, мы дополнительно воспользуемся аргументом split этой функции (подробнее см. ?summary.aov):
summary(M5, split = list(spray = list("Первые три против остальных" = 1, "ABF против CDE" = 2))) Df Sum Sq Mean Sq F value Pr(>F) spray 5 2668.8 533.8 34.702 <2e-16 *** spray: Первые три против остальных 1 93.4 93.4 6.072 0.0164 * spray: ABF против CDE 1 2558.7 2558.7 166.349 <2e-16 *** Residuals 66 1015.2 15.4 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Числа 1 и 2 в выражении, задающем значение аргумента split, указывают столбцы матрицы контрастов con.matrix, где содержатся соответствующие весовые коэффициенты.
Из полученных результатов следует, что исследованные инсектициды в целом существенно различаются по своей эффективности (см. первую строку в таблице - Pr(>F) <2e-16). Кроме того, первые три препарата в среднем значительно отличаются по эффективности действия от трех других препаратов (вторая строка: Pr(>F) = 0.0164), а препараты А, B и F в среднем отличаются от группы препаратов C, D и E (третья строка: Pr(>F) <2e-16).
---
До настоящего времени мы рассматривали только однофакторный дисперсионный анализ. Следующее сообщение будет посвящено дисперсионному анализу с двумя факторами.
Отправить комментарий