Дисперсионный анализ как частный случай общей линейной модели

В предыдущем сообщении было показано, что дисперсионный анализ (ANOVA) можно рассматривать  как линейную статистическую модель. Более того, было отмечено, что ANOVA является частным случаем т.н. общей линейной модели (ОЛМ) (General Linear Model). Понимание концепции ОЛМ очень важно для осмысленного использования lm() и других функций R, позволяющих создавать линейные модели. Поэтому стоит остановиться на ОЛМ более подробно.

Общую линейную модель можно выразить следующим образом:

Уравнение 1:

\[\mathbf{Y} = \mathbf{XB} + \mathbf{U},\]

где \(\mathbf{Y}\) - многомерная матрица предсказываемых моделью значений; \(\mathbf{X}\) - т.н. матрица модели (model matrix) или матрица плана (design matrix), которая содержит значения включенных в модель предикторов; \(\mathbf{B}\) - матрица параметров модели, которые оцениваются на основе имеющихся данных, а \(\mathbf{U}\) - матрица остатков модели.

В случае однофакторного дисперсионного анализа мы имеем дело с одномерной зависимой переменной, так что матрицы \(\mathbf{Y}\), \(\mathbf{B}\) и \(\mathbf{U}\) в приведенном уравнении (1) будут по сути представлять собой векторы из n значений каждый (где n - это объем наблюдений), тогда как размер матрицы модели \(\mathbf{X}\) будет определяться количеством уровней анализируемого фактора и объемом наблюдений. Рассмотрим этот принцип на использованном ранее примере с растениями томатов, которые выращивали при разных экспериментальных условиях - на воде (контроль; water), в среде с добавлением удобрения (Nutrient), а также в среде с добавлением удобрения и гербицида (Nutrient+24D):

# Создадим таблицу с данными:
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 (подробнее см. здесь):

tomato$trt <- relevel(tomato$trt, ref = "Water")

Графически имеющиеся данные можно представить в виде одномерной диаграммы рассеяния:

with(tomato, stripchart(weight ~ trt, xlab = "Вес, кг", ylab = "Условия"))

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

Уравнение 2:

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

где \(y_i\) - это значение веса i-го растения, \(\beta_0\) - среднее значение веса растений для базового уровня изучаемого фактора (Water), \(\beta_1\) и \(\beta_2\) - коэффициенты, отражающие разницу между средним значением веса растений в контрольной группе и средними значениями в других группах (Nutrient и Nutrient+24D соответственно), а \(\epsilon_i\) - остатки. Ключевой момент здесь состоит в том, что уровни изучаемого фактора закодированы при помощи дополнительной переменной-индикатора x (англ. dummy variable), которая может принимать только два значения - 0 или 1, в зависимости от того, к какому уровню принадлежит конкретное наблюдение \(y_i\).

Следуя концепции ОЛМ, мы можем переписать уравнение (2) в матричном виде:

Уравнение 3:

\[\begin{bmatrix}
y_1\\
y_2\\
y_3\\
y_4\\
y_5\\
y_6\\
y_7\\
y_8\\
y_9\\
y_{10}\\
y_{11}\\
y_{12}\\
y_{13}\\
y_{14}\\
y_{15}\\
y_{16}\\
y_{17}\\
y_{18}
\end{bmatrix} =
\begin{bmatrix}
1 & 0 & 0\\
1 & 0 & 0\\
1 & 0 & 0\\
1 & 0 & 0\\
1 & 0 & 0\\
1 & 0 & 0\\
1 & 1 & 0\\
1 & 1 & 0\\
1 & 1 & 0\\
1 & 1 & 0\\
1 & 1 & 0\\
1 & 1 & 0\\
1 & 0 & 1\\
1 & 0 & 1\\
1 & 0 & 1\\
1 & 0 & 1\\
1 & 0 & 1\\
1 & 0 & 1\\
\end{bmatrix} \times
\begin{bmatrix}
\beta_0 \\
\beta_1 \\
\beta_2
\end{bmatrix}\]

В уравнении (3) по левую сторону от знака "=" мы видим все измеренные значения веса растений из соответствующих экспериментальных групп (сначала идут 6 значений из контрольной группы, затем 6 значений из группы Nutrient и, наконец, 6 значений из группы Nurtient+24D). Сразу после знака "=" представлена матрица модели размером 18 x 3. Как было отмечено выше, размерность матрицы модели определяется объемом наблюдений (n = 18) и числом подлежащих оценке параметров модели (в нашем случае p = 3). Матрица модели содержит значения переменной-индикатора \(x\), которая принимает только два значения - 0 и 1, и служит для обозначения принадлежности каждого наблюдения \(y_i\) к той или иной экспериментальной группе. Матрица модели умножается на матрицу с параметрами модели, которая в рассматриваемом примере представлена вектором из 3 величин - \(\beta_0\), \(\beta_1\) и \(\beta_2\).

Согласно правилам умножения матриц, каждое значение \(y_i\) получается путем поэлементного произведения чисел из i-й строки матрицы модели на вектор, содержащий параметры модели. Например, первое значение веса из контрольной группы мы получим следующим образом:

\[y_1 = \beta_0 \times 1 + \beta_1 \times 0 + \beta_2 \times 0\ = \beta_0\]

Аналогично, первое значение веса из группы Nutrient будет рассчитано как

\[y_7 = \beta_0 \times 1 + \beta_1 \times 1 + \beta_2 \times 0\ = \beta_0 + \beta_1,\]

а первое значение веса из группы Nurtient+24D как

\[y_{13} = \beta_0 \times 1 + \beta_1 \times 0 + \beta_2 \times 1\ = \beta_0 + \beta_2\]

Очевидно, что при большом объеме наблюдений и числе параметров модели, полная запись модели по образцу уравнения (3) будет неэкономичной (представьте себе, например, ситуацию с 1000 измеренных значений зависимой переменной и 10 параметрами модели!). В связи с этим в статистике и применяют компактную матричную форму записи линейных моделей, приведенную выше в уравнении (1). Универсальность матричной формы состоит еще и в том, что матрица модели может включать значения предикторов любого типа - как количественных, так и качественных.  Так, в случае дисперсионного анализа мы имеем дело только с качественными переменными-предикторами и матрица модели содержит только значения 1 и 0, отражающие принадлежность того или иного наблюдения к конкретной группе. В случае множественной регрессии, матрица модели состоит из значений нескольких количественных предикторов. Наконец, при наличии как количественных, так и качественных предикторов мы имеем дело с т.н. ковариационным анализом (analysis of covariance, ANCOVA), и матрица модели содержит значения предикторов обоих типов. Вне зависимости от того, о каком из вариантов ОЛМ идет речь, модель всегда может быть записана в виде уравнения (1). О множественной регрессии и ковариационном анализе я расскажу в будущих сообщениях. Пока же вернемся к дисперсионному анализу, в частности к модели, представленной в уравнении (2).

Как было показано ранее, мы можем оценить параметры этой модели (\(\beta_0\), \(\beta_1\) и \(\beta_2\)) при помощи функции lm(), входящей в базовую версию R:

# Сохраним модель в виде объекта с именем М:
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

Как видим, рассчитанные параметры модели из уравнения (2) составляют \(\beta_0 = 1.683\), \(\beta_1 = 0.067\) и \(\beta_2 = -0.358\). Таким образом, модель предсказывает, что в контрольной группе средний вес растений составит 1.683 кг, в группе Nutrient он будет на 0.067 кг выше (1.683 + 0.067 = 1.750 кг), а в группе Nurtient+24D) на 0.358 кг ниже (1.683 - 0.358 = 1.325 кг), чем контрольной группе. С учетом полученных оценок параметров модели, мы можем теперь переписать уравнение (3) следующим образом:

Уравнение 4:

\[\begin{bmatrix}
y_1\\
y_2\\
y_3\\
y_4\\
y_5\\
y_6\\
y_7\\
y_8\\
y_9\\
y_{10}\\
y_{11}\\
y_{12}\\
y_{13}\\
y_{14}\\
y_{15}\\
y_{16}\\
y_{17}\\
y_{18}
\end{bmatrix} =
\begin{bmatrix}
1 & 0 & 0\\
1 & 0 & 0\\
1 & 0 & 0\\
1 & 0 & 0\\
1 & 0 & 0\\
1 & 0 & 0\\
1 & 1 & 0\\
1 & 1 & 0\\
1 & 1 & 0\\
1 & 1 & 0\\
1 & 1 & 0\\
1 & 1 & 0\\
1 & 0 & 1\\
1 & 0 & 1\\
1 & 0 & 1\\
1 & 0 & 1\\
1 & 0 & 1\\
1 & 0 & 1\\
\end{bmatrix} \times
\begin{bmatrix}
1.683\\
0.067\\
-0.358
\end{bmatrix}\]

Согласно приведенному уравнению (4), модель предсказывает, например, что первое значение веса из контрольной группы будет равно:

\[y_1 = 1.683 \times 1 + 0.067 \times 0 -0.358 \times 0 = 1.683 \]

Аналогично, первое значение веса из группы группы Nutrient будет рассчитано как

\[y_7 = 1.683 \times 1 + 0.067 \times 1 -0.358 \times 0\ = 1.683 + 0.067 = 1.750,\]
а первое значение веса из группы Nurtient+24D как

\[y_{13} = 1.683 \times 1 + 0.067 \times 0 - 0.358 \times 1\ = 1.683 - 0.358 = 1.325\]

Конечно, реальные измеренные значения веса будут в большинстве случаев несколько отличаться от предсказываемых моделью средних значений для каждой группы. Отсюда мы и получаем то, что называется остатками - это разности между наблюдаемыми и предсказанными значениями, которые не могут быть объяснены включенными в модель предикторами. Для рассматриваемого примера остатки составят:

\[ \begin{matrix} 1.5 - 1.683 = -0.183\\
1.9 - 1.683 = 0.217\\
1.3 - 1.683 = -0.383\\
1.5 - 1.683 = -0.183\\
2.4 - 1.683 = 0.717\\
1.5 - 1.683 = -0.183\\
1.5 - (1.683 + 0.067) = -0.250\\
1.2 - (1.683 + 0.067) = -0.550\\
1.2 - (1.683 + 0.067) = -0.550\\
2.1 - (1.683 + 0.067) = 0.350\\
2.9 - (1.683 + 0.067) = 1.150\\
1.6 - (1.683 + 0.067) = -0.150\\
1.9 - (1.683 - 0.358) = 0.575\\
1.6 - (1.683 - 0.358) = 0.275\\
0.8 - (1.683 - 0.358) = -0.525\\
1.15 - (1.683 - 0.358) = -0.175\\
0.9 - (1.683 - 0.358) = -0.425\\
1.6 - (1.683 - 0.358) = 0.275 \end{matrix} \]

В завершение этого сообщения стоит упомянуть о функции model.matrix(), которая позволяет просмотреть матрицу той или иной подогнанной линейной модели. Так, для созданной нами выше модели M получаем:

   (Intercept) trtNutrient trtNutrient+24D
1            1           0               0
2            1           0               0
3            1           0               0
4            1           0               0
5            1           0               0
6            1           0               0
7            1           1               0
8            1           1               0
9            1           1               0
10           1           1               0
11           1           1               0
12           1           1               0
13           1           0               1
14           1           0               1
15           1           0               1
16           1           0               1
17           1           0               1
18           1           0               1
attr(,"assign")
[1] 0 1 1
attr(,"contrasts")
attr(,"contrasts")$trt
[1] "contr.treatment"

Как видим, эта матрица в точности соответствуют матрице модели из уравнения (3). В описании атрибутов матрицы можно найти выражение "contr.treatment". О значении этого важного выражения я расскажу в следующем сообщении.