Двухфакторный дисперсионный анализ

Как следует из названия, задача рассмотренного нами ранее однофакторного дисперсионного анализа заключается в выяснении влияния какого-то одного фактора на интересующую нас количественную переменную. Однако очень редко тот или иной процесс определяется только одним фактором. Напротив - обычно  наблюдается одновременное влияние многих факторов. Задача исследователя - выявить, какие факторы оказывают существенное влияние на изучаемое явление, а какие - можно исключить из рассмотрения. Как будет показано ниже, двухфакторный дисперсионный анализ (англ. two-way analysis of variance, или two-way ANOVA) позволяет установить одновременное влияние двух факторов, а также взаимодействие между этими факторами. При наличии более двух факторов говорят о многофакторном дисперсионном анализе (англ. multifactor ANOVA; не путать с MANOVA - multivariate ANOVA!).

В качестве примера рассмотрим результаты реального эксперимента (Hand et al. 1993), в котором лабораторным крысам примерно одинакового возраста и веса в течение определенного времени давали корм с разным содержанием белка (фактор type с двумя уровнями: низкое содержание - Low, и высокое содержание - High). Кроме того, корма различались по происхождению белка (фактор source с двумя уровнями: beef - говядина, и cereal - злаки). В конце эксперимента был измерен прирост веса у крыс (weightgain) в каждой из этих групп. Таблица с данными из этого эксперимента входит в состав пакета HSAUR2, который является приложением к одной из лучших (на мой взгляд) вводных книг по R - Everitt & Hothorn (2010). Для установки этого пакет достаточно выполнить команду:

install.packages("HSAUR2")

Структура таблицы weightgain выглядит так:

library(HSAUR2)
data(weightgain)
 
str(weightgain)
'data.frame': 40 obs. of  3 variables:
 $ source    : Factor w/ 2 levels "Beef","Cereal": 1 1 1 1 1 1 1 1 1 1 ...
 $ type      : Factor w/ 2 levels "High","Low": 2 2 2 2 2 2 2 2 2 2 ...
 $ weightgain: int  90 76 90 64 86 51 72 90 95 78 ...

Разведочный анализ данных

Перед выполнением любого статистического анализа полезно "рассмотреть" данные на графике (подробнее об использованном пакете ggplot2 см. здесь):

library(ggplot2)
 
ggplot(data = weightgain, aes(x = type, y = weightgain)) + 
       geom_boxplot(aes(fill = source))

Также стоит ознакомиться со сводными описательными статистиками (например, используя возможности пакета doBy):

require(doBy)
 
summaryBy(weightgain ~ type + source, data = weightgain,
          FUN = c(mean, sd, length))
  type source weightgain.mean weightgain.sd weightgain.length
1 High   Beef           100.0      15.13642                10
2 High Cereal            85.9      15.02184                10
3  Low   Beef            79.2      13.88684                10
4  Low Cereal            83.9      15.70881                10

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

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

plot.design(weightgain)

Из полученного графика видно, что наибольшая разница в средних приростах веса крыс связана с уровнем содержания белка в корме, тогда как эффект источника происхождения белка выражен в меньшей степени.

Рассматриваемый эксперимент мы можем отнести к т.н. полнофакторному эксперименту (англ. full factorial experiment), поскольку в нем реализуются все возможные сочетания имеющихся уровней факторов. Значительное преимущество такого дизайна эксперимента заключается в том, что он позволяет выяснить наличие взаимодействия между изучаемыми факторами. В рамках дисперсионного анализа, под "взаимодействием" (англ. interaction) понимают такую ситуацию, когда переменная-отклик ведет себя по разному при разных сочетаниях изучаемых факторов. Понять эту концепцию поможет "график взаимодействий" (interaction plot), который в R можно построить при помощи базовой функции interaction.plot():

with(weightgain, interaction.plot(x.factor = type,
                                  trace.factor = source,
                                  response = weightgain))

Из приведенного рисунка видно, что при высоком содержании белка в корме, прирост веса крыс в среднем также высок, но при условии, что этот белок имеет животное происхождение. Если же содержание белка низкое, то ситуация меняется на противоположную - прирост оказывается несколько выше (хотя и не намного) в группе крыс, получавших корм растительного происхождения.


Выполнение дисперсионного анализа при помощи функции aov()

Помимо того, что рассматриваемый эксперимент является полнофакторным, во всех четырех группах имеется также одинаковое число крыс (по 10 в каждой), т.е. мы имеем дело со сбалансированным набором данных. Как было показано ранее, для анализа сбалансированных наборов данных мы можем применить классический способ разложения общей дисперсии в данных на отдельные составляющие, реализованный в функции aov():

M1 <- aov(weightgain ~ source + type + source:type, data = weightgain)
 
summary(M1)
            Df Sum Sq Mean Sq F value Pr(>F)  
source       1    221   220.9   0.988 0.3269  
type         1   1300  1299.6   5.812 0.0211 *
source:type  1    884   883.6   3.952 0.0545 .
Residuals   36   8049   223.6                 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Подробное объяснение того, как интерпретировать полученную таблицу дисперсионного анализа, было приведено ранее. В целом, можно сделать вывод об отсутствии статистически значимой связи между приростом веса крыс и источником белка в корме (P = 0.3269), тогда как влияние уровня содержания белка в корме оказалось значимым (Р = 0.0211). Взаимодействие между источником происхождения белка и уровнем его содержания в корме незначимо (P = 0.0545), что, в принципе, согласуется с результатом анализа приведенного выше графика взаимодействий.

Обратите внимание на то, как в формуле, поданной на функцию aov(), было задано взаимодействие между двумя факторами: сначала были приведены два главных фактора, разделенные знаком "+", а затем к ним добавлено выражение "source:type". Это стандартный синтаксис для такого рода анализа в R. Однако приведенную формулу можно было бы также сократить до weightgain ~ source*type - результат оказался бы идентичным. Подробнее о синтаксисе формул в R можно узнать, например, здесь (англ. яз.).


Выполнение дисперсионного анализа при помощи функции lm()

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

Применим функцию lm() в отношении данных по приросту веса у крыс:

M2 <- lm(weightgain ~ type*source, data = weightgain)
 
summary(M2)
Call:
lm(formula = weightgain ~ type * source, data = weightgain)
 
Residuals:
   Min     1Q Median     3Q    Max 
-29.90  -9.90   2.05  10.85  25.10 
 
Coefficients:
                     Estimate Std. Error t value Pr(>|t|)    
(Intercept)           100.000      4.729  21.148  < 2e-16 ***
typeLow               -20.800      6.687  -3.110  0.00364 ** 
sourceCereal          -14.100      6.687  -2.109  0.04201 *  
typeLow:sourceCereal   18.800      9.457   1.988  0.05447 .  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
 
Residual standard error: 14.95 on 36 degrees of freedom
Multiple R-squared:   0.23, Adjusted R-squared:  0.1658 
F-statistic: 3.584 on 3 and 36 DF,  p-value: 0.02297

В соответствии с рассмотренными ранее принципами, полученные параметры модели мы интерпретируем следующим образом. В первой строке таблицы c параметрами модели (Intercept) представлена информация, относящаяся к среднему значению прироста веса в группе крыс, которым давали корм с высоким содержанием белка (High) животного происхождения (Beef) (для простоты обозначим эту базовую группу как "High - Beef"). Видим, что средний прирост веса в этой группе составил 100 г, и что этот прирост значимо отличается от 0 (стандартная ошибка = 4.729 г, значение t-критерия Стьюдента = 21.148, P-значение << 0.001). Вторая строка (typeLow) содержит оценку величины эффекта, связанного с низким содержанием в корме белка животного происхождения: в этой группе крыс прирост веса оказался существенно ниже (на 20.8 г, Р = 0.00364), чем в группе "High - Beef" (100 - 20.8 = 79.2 г). Из третьей строки таблицы (sourceCereal) мы узнаём, что у крыс, которым давали корм с высоким содержанием белка растительного поисхождение прирост веса был статистически значимо ниже (на 14.1 г, Р = 0.04201), чем в группе "High - Beef" (100 - 14.1 = 85.9 г). Наконец, в четвертой строке мы находим оценку эффекта взаимодействия между факторами "type" и "source" - 18.8 г. Сложив первые три коэффициента с эффектом взаимодействия, получим среднее значение прироста веса крыс в группе, которым давали корм с низким содержанием белка растительного происхождения: 100 - 20.8 - 14.1 + 18.8 = 83.9 (см. также приведенную ранее таблицу со средними значениями по каждой группе). Как видим, этот средний прирост лишь незначительно (Р = 0.05447) превысил таковой в группе крыс, которым давали корм с низким содержанием белка растительного происхождения (83.9 vs. 79.2). Как видим, полученные результаты хорошо согласуются с впечатлениями, которые мы получили в ходе расчета описательных статистик по каждой группе, а также при выполнении графического разведочного анализа данных (см. выше).

Как было показано ранее, результаты дисперсионного анализа, выполненного при помощи функции lm(), можно представить также и в виде классической ANOVA-таблицы. Для этого соответствующую модель необходимо подать на функцию anova() (не путать с aov()!):

anova(M2)
Analysis of Variance Table
 
Response: weightgain
            Df Sum Sq Mean Sq F value  Pr(>F)  
type         1 1299.6 1299.60  5.8123 0.02114 *
source       1  220.9  220.90  0.9879 0.32688  
type:source  1  883.6  883.60  3.9518 0.05447 .
Residuals   36 8049.4  223.59                  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Представленные в этой таблице результаты идентичны таковым в таблице, полученной ранее при помощи функции aov().


О порядке перечисления факторов в формуле модели

В ходе спецификации модели M2 мы сначала указали фактор "type", а затем "source" (см. выше: weightgain ~ type*source). Что произойдет, если этот порядок изменить?

M3 <- lm(weightgain ~ source*type, data = weightgain)
 
anova(M3)
Analysis of Variance Table
 
Response: weightgain
            Df Sum Sq Mean Sq F value  Pr(>F)  
source       1  220.9  220.90  0.9879 0.32688  
type         1 1299.6 1299.60  5.8123 0.02114 *
source:type  1  883.6  883.60  3.9518 0.05447 .
Residuals   36 8049.4  223.59                  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Результаты идентичны тем, что были получены для модели M2. Единственная разница состоит лишь в том, что строки, содержащие информацию по факторам "type" и "source" поменялись местами.  Однако, как было отмечено выше, такое полное совпадение результатов будет наблюдаться только при работе со сбалансированными наборами данных. Для демонстрации этого утверждения изменим исходные данные путем удаления, например, первых 6 наблюдений и последних 7 наблюдений:

set.seed(123)
weightgain2 <- weightgain[-c(1:6, 34:40), ]
 
# Модели с разным порядком указания предикторов:
M4 <- lm(weightgain ~ type*source, data = weightgain2)
M5 <- lm(weightgain ~ source*type, data = weightgain2)
 
# ANOVA-таблица для модели М4
anova(M4)
Analysis of Variance Table
 
Response: weightgain
            Df Sum Sq Mean Sq F value  Pr(>F)  
type         1  758.0  758.02  3.1655 0.08843 .
source       1  584.8  584.76  2.4419 0.13179  
type:source  1  744.5  744.54  3.1092 0.09113 .
Residuals   23 5507.6  239.46                  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
 
# ANOVA-таблица для модели М5
anova(M5)
Analysis of Variance Table
 
Response: weightgain
            Df Sum Sq Mean Sq F value  Pr(>F)  
source       1 1188.8 1188.83  4.9645 0.03594 *
type         1  153.9  153.95  0.6429 0.43087  
source:type  1  744.5  744.54  3.1092 0.09113 .
Residuals   23 5507.6  239.46                  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Как видим, выводы, сделанные на основании этих двух новых моделей оказываются совершенно разными: ни один из параметров модели М4 не является статистически значимым, тогда как модель M5 указывает на существование существенного эффекта источника происхождения белка. Такая разница обусловлена алгоритмом, по которому происходит декомпозиция общей дисперсии при анализе несбалансированных наборов данных. Как же быть? Как из последних двух моделей отражает действительность?

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

summary(M4)
 
Call:
lm(formula = weightgain ~ type * source, data = weightgain2)
 
Residuals:
   Min     1Q Median     3Q    Max 
-27.00 -10.82   2.00  11.18  23.10 
 
Coefficients:
                     Estimate Std. Error t value Pr(>|t|)    
(Intercept)           100.000      4.893  20.435 3.02e-16 ***
typeLow               -16.250      9.155  -1.775   0.0891 .  
sourceCereal          -24.000     10.187  -2.356   0.0274 *  
typeLow:sourceCereal   24.150     13.696   1.763   0.0911 .  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
 
Residual standard error: 15.47 on 23 degrees of freedom
Multiple R-squared:  0.2748, Adjusted R-squared:  0.1802 
F-statistic: 2.906 on 3 and 23 DF,  p-value: 0.05643
 
 
summary(M5)
 
Call:
lm(formula = weightgain ~ source * type, data = weightgain2)
 
Residuals:
   Min     1Q Median     3Q    Max 
-27.00 -10.82   2.00  11.18  23.10 
 
Coefficients:
                     Estimate Std. Error t value Pr(>|t|)    
(Intercept)           100.000      4.893  20.435 3.02e-16 ***
sourceCereal          -24.000     10.187  -2.356   0.0274 *  
typeLow               -16.250      9.155  -1.775   0.0891 .  
sourceCereal:typeLow   24.150     13.696   1.763   0.0911 .  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
 
Residual standard error: 15.47 on 23 degrees of freedom
Multiple R-squared:  0.2748, Adjusted R-squared:  0.1802 
F-statistic: 2.906 on 3 and 23 DF,  p-value: 0.05643

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

Поскольку источник происхождения теперь имеет больший эффект на прирост веса крыс, чем уровень содержания белка, предпочтение следует отдать модели М5.


Многофакторный дисперсионный анализ

Конечно, ничто не запрещает нам включать в анализ больше, чем два фактора. Тщательно спланированный и выполненный многофакторный дисперсионный анализ может оказаться очень полезным для понимания изучаемого явления. Однако следует помнить о возникающих при этом затруднениях. Уже при наличии трех основных факторов (X, Y и Z) в модель войдут четыре дополнительных параметра, отражающих взаимодействия между этими факторами (XY, XZ, YZ, XYZ). Чем больше параметров мы оцениваем, тем больше степеней свободы нам нужно. Другими словами, требуется собрать больше данных, что не всегда возможно. Кроме того, при проверке статистических гипотез в отношении большого числа параметров возрастает вероятность сделать неверное заключение о значимости как минимум некоторых из них. Но даже если мы имеем "достаточно" данных для надежной оценки большого числа параметров, возникает также проблема с интерпретацией этих параметров - в частности, с интерпретацией взаимодействий высокого порядка (>2-3). К счастью, однако, в большинстве случаев взаимодействия высокого порядка оказываются незначимыми, что дает нам возможность исключить их из модели. Выбор оптимальных моделей - это отдельная большая тема, к которой я обязательно вернусь в будущем.