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


Суть "классического" дисперсионного анализа мы рассмотрели на примере эксперимента (Maindonald & Braun 2010), в котором томаты выращивали при трех разных экспериментальных условиях (на воде, в среде с добавлением удобрения, а также в среде с добавлением удобрения и гербицида). В конце эксперимента был измерен вес растений из этих трех групп. Мы выяснили, что каждое измеренное значение веса растений можно разложить на следующие составляющие:

Уравнение 1:  \[x_{ij} = \bar{X} + (\bar{x_i}- \bar{X}) + (x_{ij} - \bar{x_i})\]

где \((\bar{x_i}- \bar{X})\) - отклонения групповых средних от общего среднего значения (т.е. вычисленного по всем значениям веса растений), а \((x_{ij} - \bar{x_i})\) - отклонения отдельных наблюдений от среднего значения группы, к которой они принадлежат. Используя эти отклонения, а также учитывая число анализируемых групп и число наблюдений в каждой группе, мы далее рассчитали меж- и внутригрупповую дисперсии и сравнили их при помощи F-критерия Фишера. Исходя из величины полученного критерия, был сделан вывод о вероятности справедливости нулевой гипотезы об отсутствии разницы между групповыми средними значениями веса растений.

Приведенное выше уравнение (1) для однофакторного дисперсионного анализа мы можем записать в более простом виде:

Уравнение 2:

\[y_{ij} = \mu + \alpha_{i} + \epsilon_{ij}\]

где \(y_{ij}\) - вес растения \(j\) из экспериментальной группы \(i\), \(\mu\) - среднее значение веса всех задействованных в эксперименте растений, \(\alpha_{i}\) - разность между \(\mu\) и средним значением растений в группе \(i\), a \(\epsilon_{ij}\) - это остатки, отражающие отклонения отдельных наблюдений в группе \(i\) от среднего значения этой группы. Предполагается, что остатки имеют нормальное распределение со средним значением 0 и стандартным отклонением \(\sigma\), причем это стандартное отклонение одинаково во всех группах.

Как можно заметить, уравнение (2) почти ничем не отличается от уравнения простой регрессии, которая, выражаясь в терминах статистики, является линейной моделью (прилагательное "линейная" подчеркивает аддитивный характер эффектов, оказываемых предикторами (переменные по правую сторону от знака "равно") на переменную-отклик):

Уравнение 3:

\[y_{i} = \beta_0 + \beta_1 x_{i} + \epsilon_{i}\]

где \(y_i\) - это i-е значение зависимой количественной переменной \(y\); \(x_i\) - значение  количественного предиктора \(x\), соответствующее i-тому значению \(y\); \(\beta_0\) - значение, которое принимает \(y\) при \(x = 0\) (геометрически, это точка, в которой линия регрессии пересекает ось Y; отсюда англ. термин для этого коэффициента - intercept, что значит "пересечение"); \(\beta_1\) - коэффициент регрессии, который показывает, насколько изменяется \(y\) при изменении \(x\) на одну единицу; \(\epsilon_{i}\) - это остатки, т.е. разница между наблюдаемыми значениями \(y_i\) и средними значениями, предсказанными моделью для каждого \(x_i\).

Наиболее важные отличия между уравнениями (2) и (3) заключаются в следующем:
  • В обоих случаях зависимая переменная \(y\) (=переменная-отклик) является количественной. Однако в дисперсионном анализе значения \(y\) сгруппированы в соответствии с \(i\) уровнями изучаемого фактора - отсюда индекс \(i\) в обозначении \(y_{ij}\), указывающий, к какой группе принадлежит каждое наблюдение.
  • В ANOVA независимая переменная (=предиктор) является качественной, тогда как в регрессионной модели мы имеем дело с количественным предиктором. Степень взаимосвязи между предиктором и независимой переменной в регрессионной модели определяется коэффициентом регрессии (\(\beta_1\) в уравнении (3)). В дисперсионном же анализе размер эффекта изучаемого фактора (т.е., степень отклонения групповой средней от общего среднего значения) выражен при помощи \(\alpha_{i}\).
К сожалению, с моделью в уравнении (2) имеется одна загвоздка - не все ее параметры  идентифицируемы. Это значит, что существует бесконечное множество возможных (неуникальных) значений параметров модели. Например, если мы добавим к \(\mu\) 2.5 (или любую другую константу) и соответственно вычтем 2.5 из каждого \(\alpha_i\), то получим точно те же значения \(y_{ij}\). Однако проблему неидентифицируемости можно легко обойти, наложив определенные ограничения на параметры модели.

При работе с категориальными предикторами (факторами) возможны несколько разных подходов к наложению подобных ограничений (см., например, Faraway 2009, с. 191-192). Во многих статистических программах, включая R, величину эффекта для первого уровня анализируемого фактора ("базового" уровня; англ. reference level или baseline level) по умолчанию принимают равной 0: \(\alpha_1 = 0\). Тогда \(\mu\) будет представлять собой среднее значение зависимой переменной для базового уровня фактора (например, в контрольной группе какого-либо эксперимента), тогда как \(\alpha_i\) для \(i \neq 1\) будут отражать разницу между средним значением базового уровня и средними значениями остальных уровней. Используя этот подход, мы можем переписать уравнение (2) в виде регрессионного уравнения следующим образом:

Уравнение 4:

\[y_{i} = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \dots + \beta_{k-1} x_{k-1} + \epsilon_{i}\]

где \(\beta_0\) - среднее значение зависимой переменной для базового уровня изучаемого фактора (например, в контрольной группе); \(\beta_1 \dots \beta_{k-1}\) - коэффициенты, отражающие разницу между средним значением базового уровня и средними значениями остальных уровней; \(\epsilon_i\) - остатки. Ключевой момент для понимания здесь состоит в том, что уровни изучаемого фактора (всего имеется k таких уровней) закодированы при помощи дополнительной переменной-индикатора \(x\) (англ.  dummy variable, дословно "переменная-пустышка"), которая может принимать только два значения - 0 или 1, в зависимости от того, к какому уровню принадлежит конкретное наблюдение \(y_i\). Например, если наблюдение \(y_3\) принадлежит к уровню \(x_1\), тогда \(x_1\) в уравнении (4) примет значение 1, а \(x_2 \dots x_{k-1}\) все будут равны 0, и следовательно для наблюдения \(y_3\) регрессионное уравнение примет вид:

\[y_3 = \beta_0 + \beta_1 \times 1 + \epsilon_3 = \beta_0 + \beta_1 + \epsilon_3\]

Уравнение (4) как раз и описывает ту линейную модель, которая рассчитывается в R при выполнении однофакторного дисперсионного анализа с помощью функции lm(). Посмотрим, как это всё работает, вернувшись к примеру с растениями томатов, которые выращивали при разных экспериментальных условиях.

# Создадим таблицу с данными:
tomato <-
  data.frame(weight=
               c(1.5, 1.9, 1.3, 1.5, 2.4, 1.5, # water
                 1.5, 1.2, 1.2, 2.1, 2.9, 1.6, # nutrient
                 1.9, 1.6, 0.8, 1.15, 0.9, 1.6), # nutrient+24D
             trt = rep(c("Water", "Nutrient", "Nutrient+24D"),
                       c(6, 6, 6)))

В зависимости от стоящей задачи, исследователь волен выбрать в качестве базового любой  уровень изучаемого фактора. В рассматриваемом примере имеет смысл проводить сравнения с уровнем Water (растения, выращенные на воде). Если базовый уровень не задан исследователем однозначно, R автоматически расположит названия всех уровней в алфавитном порядке и выберет в качестве базового первый уровень из этого списка (например, в нашем случае это был бы уровень Nutrient). Задать базовый уровень фактора в R позволяет функция relevel():

# Автоматически R выберет Nutrient в качестве базового уровня:
levels(tomato$trt)
[1] "Nutrient"     "Nutrient+24D" "Water"
 
# Изменим базовый уровень на Water:
tomato$trt <- relevel(tomato$trt, ref = "Water")
 
levels(tomato$trt)
[1] "Water"        "Nutrient"     "Nutrient+24D"

Теперь рассчитаем коэффициенты линейной модели для рассматриваемого эксперимента:

# Сохраним модель в виде объекта с именем М:
M <- lm(weight ~ trt, data = tomato)
 
# Просмотрим результаты вычислений:
summary(M)
 
Call:
lm(formula = weight ~ trt, data = tomato)
 
Residuals:
    Min      1Q  Median      3Q     Max 
-0.5500 -0.3500 -0.1792  0.2750  1.1500 
 
Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)      1.68333    0.20849   8.074 7.69e-07 ***
trtNutrient      0.06667    0.29485   0.226    0.824    
trtNutrient+24D -0.35833    0.29485  -1.215    0.243    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
 
Residual standard error: 0.5107 on 15 degrees of freedom
Multiple R-squared: 0.1381, Adjusted R-squared: 0.02321 
F-statistic: 1.202 on 2 and 15 DF,  p-value: 0.328



Разберемся с основными элементами приведенной таблицы с результатами анализа. Прежде всего, программа напоминает нам, к какой модели относятся выводимые результаты (это бывает очень удобно при одновременном построении большого количества моделей): Call: lm(formula = weight ~ trt, data = tomato). Далее следует сводная информация о распределении остатков (Residuals: ... ; на эту тему мы поговорим позднее).

Самая интересная часть таблицы с результатами анализа представлена под заголовком Coefficients. Здесь приведены оценки коэффициентов построенной нами модели

\[y_{i} = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \epsilon_{i}\]

В первой строке (Intercept) мы видим информацию, относящуюся к коэффициенту \(\beta_0\). В столбце Estimate ("оценка") представлено рассчитанное значение этого коэффициента (1.683), которое равно среднему весу растений в группе Water (установленный нами базовый уровень изучаемого фактора). Во второй строке (trtNutrient) мы находим коэффициент \(\beta_1\), который отражает разность между базовым уровнем и средним значением в группе Nutrient (1.683 + 0.067 = 1.750). Наконец, в третьей строке (trtNutrient+24D) мы видим коэффициент \(\beta_2\), который отражает разность между базовым уровнем и средним значением в группе Nutrient+24D (1.683 - 0.358 = 1.325). Верность описанной интерпретации полученных коэффициентов мы можем легко проверить, рассчитав средние значения для каждой экспериментальной группы:

tapply(weight, trt, mean)
       Water     Nutrient Nutrient+24D 
    1.683333     1.750000     1.325000

В столбцах Std. Error, t value и Pr(>|t|) представлены стандартные ошибки рассчитанных коэффициентов, значения t-критерия Cтьюдента (используется для проверки гипотезы о равенстве коэффициента 0) и Р-значения для каждого t-критерия соответственно. Как видим, ни один из исследованных уровней фактора "условия выращивания" не отличается от базового уровня по среднему значению веса растений - в обоих случаях Р-значения регрессионных коэффициентов оказались >0.05, что указывает на то, что наблюдаемое их отличие от нуля с большой долей вероятности является случайным (другими словами, исследованные условия выращивания в среднем не оказали никакого влияния на вес растений по сравнению с растениями, выращенными на воде).

Далее в таблице с результатами анализа мы видим:
  • Residual standard error - оценка стандартного отклонения остатков. Вспомним: предполагается, что остатки имеют нормальное распределение со средним значением 0 и стандартным отклонением \(\sigma\). В данной строке результатов анализа как раз и приводится оценка значения \(\sigma\).
  • Multiple R-squared и Adjusted R-squared: коэффициент детерминации и коэффициент детерминации с поправкой на число параметров модели соответственно.
  • F-statistic: значение критерия Фишера, при помощи которого проверяется нулевая гипотеза о том, что все коэффициенты модели (в нашем случае \(\beta_1\) и  \(\beta_2\)) равны 0.
Значение F-критерия для линейных моделей рассчитывается как отношение дисперсии в данных, "объясненной" параметрами модели, к дисперсии остатков (т.е., части общей дисперсии, которую модель "не объясняет"). Это легко легко увидеть, подав объект M на функцию anova() (не путать с aov()):

anova(M)
Analysis of Variance Table
 
Response: weight
          Df Sum Sq Mean Sq F value Pr(>F)
trt        2 0.6269 0.31347  1.2019  0.328
Residuals 15 3.9121 0.26081

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

summary(aov(weight ~ trt, data = tomato))
            Df Sum Sq Mean Sq F value Pr(>F)
trt          2  0.627  0.3135   1.202  0.328
Residuals   15  3.912  0.2608

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

Тесное сходство между моделью дисперсионного анализа и регрессионной моделью неслучайно, поскольку обе являются частными случаями общей линейной модели (англ. General Linear Model(!) не путать с обобщенными линейными моделями - Generalized Linear Models). Об этом классе статистических моделей, а также о понятии контрастов в дисперсионном анализе, я расскажу в следующем сообщении.


Литература:
Faraway JJ (2009) Linear models with R. Chapman & Hall/CRC


6 Комментарии

Andrey ILIN написал(а)…
Спасибо!
Unknown написал(а)…
Спасибо, очень помогли разобраться!
Анонимный написал(а)…
Сергей, скажите, а с чем сравнивается базовый уровень? Почему на против него тоже стоит оценка статистической значимости и Т-критерия Стьюдента?
Sergey Mastitsky написал(а)…
Проверяется нулевая гипотеза о равенстве соответствующего генерального среднего значения нулю.
Феофан написал(а)…
Феофан
а у меня при повторении всех ваших расчетов не срабатывает
tapply(weight,trt,mean)
Unknown написал(а)…
здравствуйте, объясните пожалуйста, что означает таблица Residuals...
Новые Старые