Про контрасты в R, которые не совсем контрасты

С контрастами сталкивались, наверное, все, кто имеет дело с линейными моделями. Тема достаточно подробно раскрыта в книге Статистический анализ и визуализация данных с помощью R (pdf; сообщение в блоге), а также в книге Р.Кабакова "R в действии" и во многих других источниках.

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


Если коротко, то "контрасты" в R заданы в соответствии со схемой записи линейной модели при помощи фиктивных переменных, принимающих значения 0, 1 или 0, 1, -1.

Разберем простой пример из книги Logan M. "Biostatistical Design and Analysis Using R".

Пусть есть такие данные:

> Y <- c(2, 3, 4, 6, 7, 8, 10, 11, 12)
> A <- gl(3, 3, 9, lab = c("G1", "G2", "G3"))

Зададим контрасты, которые R использует по умолчанию, т.е. контрасты комбинаций условий (или комбинаций воздействий, также известные как простые контрасты), которые сравнивают с первым (базовым) уровнем все остальные уровни фактора:

> contrasts(A) <- contr.treatment
> contrasts(A)
      2 3
G1 0 0
G2 1 0
G3 0 1

Можно заметить, что представленная матрица "контрастов" как таковых контрастов не содержит. Алгоритм получения данной матрицы следующий.

1. Записать оцениваемые коэффициенты линейной модели с контрастами в виде суммы средних значений сравниваемых уровней. Для фактора с тремя уровнями и простых контрастов с базовым первым уровнем:

\[\beta_0=1\mu_1+0\mu_2+0\mu_3\] \[\beta_1=-1\mu_1+1\mu_2+0\mu_3\] \[\beta_2=-1\mu_1+0\mu_2+1\mu_3\]

В матричной форме:

> mat <- rbind(c(1, 0, 0), c(-1, 1, 0), c(-1, 0, 1))
> mat
      [,1] [,2] [,3]
[1,]    1    0    0
[2,]   -1    1    0
[3,]   -1    0    1

2. Найти обратную матрицу:

> solve(mat)
     [,1] [,2] [,3]
[1,]    1    0    0
[2,]    1    1    0
[3,]    1    0    1

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

Об этом важно помнить при использовании пользовательских контрастов. Следующий пример из вышеназванной книги (пример неправильного способа задать контрасты). Сравниваем среднее значение группы G1 с общим средним значением групп G2 и G3:

> contrasts(A) <- cbind(c(1, -0.5, -0.5))
> contrasts(A)
   [,1]          [,2]
G1  1.0 -6.410345e-17
G2 -0.5 -7.071068e-01
G3 -0.5  7.071068e-01
> l <- lm(Y~A)
> summary(l)$coefficients
             Estimate Std. Error   t value     Pr(>|t|)
(Intercept)  7.000000  0.3333333 21.000000 7.595904e-07
A1          -4.000000  0.4714045 -8.485281 1.465426e-04
A2           2.828427  0.5773503  4.898979 2.713682e-03

Свободный член, равный 7, соответствует общему среднему всех трех групп. Среднее первой группы равно 3, второй - 7, третьей - 11; общее среднее второй и третьей группы - 9. Тогда второй коэффициент должен быть равен (3-7)-(9-7)=3-9=-6, а не -4, как у нас получилось.
Правильный вариант:

> mat <- rbind(c(1, 1, 1), c(1, -0.5, -0.5), c(0.5, 0.5, -1))
> mat
     [,1] [,2] [,3]
[1,]  1.0  1.0  1.0
[2,]  1.0 -0.5 -0.5
[3,]  0.5  0.5 -1.0
> mat <- solve(mat)
> mat
          [,1]       [,2]          [,3]
[1,] 0.3333333  0.6666667 -3.702550e-17
[2,] 0.3333333 -0.6666667  6.666667e-01
[3,] 0.3333333  0.0000000 -6.666667e-01
> contrasts(A) <- mat[, -1]
> l <- lm(Y~A)
> summary(l)$coefficients
            Estimate Std. Error   t value     Pr(>|t|)
(Intercept)        7  0.3333333 21.000000 7.595904e-07
A1                -6  0.7071068 -8.485281 1.465426e-04
A2                -6  0.7071068 -8.485281 1.465426e-04

Здесь второй заданный контраст сравнивает среднее третьей группы с общий средним первых двух. Он, как и первый, равен -6.

Написано по мотивам http://sites.stat.psu.edu/~jls/stat512/lectures/ContrastsInR.pdf. Возможно, кому-то в первоисточнике будет понятнее.